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Foreword 



Most of the information available in modern societies is visual. We have huge amounts of images 
and videos obtained by a wide variety of sensors: scanners, video cameras, medical imaging 
systems (e.g., CT, MRI, ultrasound, SPECT, PET), radar, sonar or satellite systems. 

Each type of images has its own properties and raises new challenges since it measures a small 
part of the world in its own way. For example different image modalities (MRI, SPECT) reveal 
complementary aspects of the patient (anatomy and function) useful for diagnosis. Satellite 
images convey information about the land cover. Despite this diversity there are common 
goals. The common goal is to understand the information conveyed by each type of image and 
to devise methods to extract this information and fuse information from multiple images. 

The areas of Image Processing and Computer Vision deal with a wide range of problems 
which can be grouped in terms of two main questions: 

Problem 1: How can we modify images in order to enhance their properties? 

Problem 2: How can we extract information from images in an automatic way? 

In the first case, the outputs are images. Given one or several input images, we wish to 
modify them in order to enhance information (e.g., remove noise, correct a geometric distortion, 
compare satellite images obtained at different years). In the second problem, the output is a 
description of the image content (e.g., roads, people, human organs, 3D shape of a scene). 

Problem 1 is the key problem in Image Processing and Problem 2 is the key problem in 
Image Analysis and Vision. The difference between Image Analysis and Vision is not always 
clear. In Image Analysis we are interested in describing the image content (e.g., detection of 
faces in images) while in Computer Vision we wish to extract 3D information about the real 
world (e.g., estimation of 3D shape and motion) and control it (e.g., robot guidance). 

A few examples are listed below. 
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Problems: 

• identify people from biometric signals (iris, fingerprint, signature, face image) 

• detect people faces in images 

• monitor urban zones and detect abnormal behaviors (surveillance) 

• align medical images and estimate the boundary of organs 

• detect and classify diseases from medical images 

• classify land cover using satellite images 

• improve the quality of degraded images 



These problems are very different and they are solved by different methods. However, all 
methods share common steps: i) development of a model for the information to be retrieved, 
ii) fitting the model to the image or group of images and iii) evaluation of the model. These 
three steps can be performed in different ways as we shall see during this course. 

There is not a unified theory for Vision or for the human perception as we have in the field of 
Mechanics with the Lagrangean framework or in Electromagnetism with the Maxwell equations. 
The best idea we can provide about this field is that we have many challenges (problems) and 
a wide variety of mathematical methods (tools) . When we want to solve a new problem we can 
probably use several methods we know. Experience allows us to guess which ones will work 
best and which kind of changes will be needed. 

Image processing and computer vision are strongly connected to several disciplines e.g., 
statistical modeling, optimization, pattern recognition and tracking. However, it is fair to say 
that many of the techniques used in image analysis and vision can be classified as statistical 
learning methods. 

This course tries to identify some of the key methods in Computer Vision and Image Process- 
ing and provides a discussion of a few challenging problems hoping that this approach will 
provide a useful introduction to the field. 
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Sources 

There are excellent textbooks in specific areas of Image Processing and Computer Vision e.g., 
3D vision [1, 2], wavelets [3], image reconstruction [4]). 

It is more difficult to suggest a single textbook with broad scope. The area is so vast 
that it is difficult to select the topics and to identify which are the most important. It all 
depends on the focus. There is also a difficult trade off between a large scope and depth. 
Some classic books which cover many aspects of image processing and computer vision are 
[6, 7, 8, 9, 10]. An good alternative for a broad introduction to this field is given by the Handbook 
of Image and Video Processing [11] written by several top researchers in this area. Electronic 
versions of several books are available in the CVOnline site (http://homepages.inf.ed.ac.uk 
/rbf/CVonline/books.htm). 

Despite the effort to write good textbooks most of the information is available on scientific 
journals. A few examples are: 

• Transactions Pattern Analysis and Machine Vision 

• Transactions on Image Processing 

• Transactions Medical Imaging 

• International Journal of Computer Vision 

• Image and Vision Computing 



The tour is about to begin. I hope you like it. 
Jorge S. Marques 
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Part I 
Basics 



Chapter 1 



Images 



1.1 What is an image? 

Images are signals which represent the evolution of a quantity in the plane or in space. 

Let us consider for example a gray level image. Each image point x is associated to a real 
number (intensity) /(x) which measures the amount of radiation received at the vicinity of x. 
The set of points for which the signal is defined in called the domain V. A gray level image is 
a signal which associates a real value to each point xGP. 

Figure 1.1 shows a well known image called boat (left). The image can be visualized as an 
intensity map which associates an elevation to each point of the domain. The intensity map 
for a small block of the image boat is shown in figure 1.1 (center). Try to guess which region 
is represented without looking at the block on the top. It is not simple. We are not used to 
interpret such intensity surfaces. The intensity surface has the same information as the original 
image. However we are not able to analyse its content most of the time since our visual system 
is not trained for that task. 

There are several types of images. Typical examples are color images. The color is a 
human perception, a response of the human brain to the radiation received at the retina cells. 
Fortunately, few information is extracted from electromagnetic spectrum by the retina and only 
three coefficients (e.g., the intensities of three primary colors RGB) are needed to specify the 
color perception. Color images are represented in a computer by the intensities of R (red), G 
(green) and B (blue) i.e. by a vector with three components. Therefore, a color image is a signal 
which associates each point of the domain with a vector belonging to M 3 , i.e., / : V — > M 3 . 
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Figure 1.1: Image boat, intensity surface (block) and intensity profile along the main diagonal 

From a mathematical point of view, an image is a map / :V — ^ 1Z where the domain V is 
a subset of M 2 or R 3 . 



Examples 

Two examples are the images defined in the interval of M 2 e.g. in a domain V — [0, l] 2 or 
in a grid of points V = {0, 1, . . . , M - 1} x {0, 1, . . . , N - 1} 



When the domain is M 2 or an interval of M 2 the image is denoted as a continuous image. 
When the domain is Z 2 or an interval of Z 2 the image is called a discrete image. 

This raises a question: is it possible to store a continuous image since it has an infinite 
number of points? we will address this question later. 

Range 

The range of an image depends on the sensor used during the acquisition process and the 
transformations applied to the image. Examples: 

• real images are produced by many sensors: CCD cameras, computed tomography, ul- 
trasound sensors, etc. In this case the co-domain is 1Z = R. 

• vector images map each point of the domain to a vector of (1Z = W 1 ). The simplest case 
are the images produced by color cameras which associate a vector with three components 



11 



(RGB) to each point. Another example are satellite images which measure the intensity 
of the radiation in several frequency bands from the visible to the infra-red. 

• complex images are produced by magnetic resonance imaging systems (MRI) and in 
many image processing operations (e.g., Fourier transform). They map each point of the 
domain to a complex number. 

• label images associate a label to each point. The label can be represented by an integer 
number belonging to a set of admissible labels 1Z — {0,1,..., L — 1}. The labels may 
identify different objects in the image (p. ex., vehicles, people, obstacles). 



1.2 Discrete Images 



Discrete images are images whose domain is Z 2 or a subset of Z 2 . Most sensors produce finite 
discrete images with M x N elements. The elements are usually associated to a rectangular 
grid and can be stored in a matrix 



1 = 



7(0,0) 
7(1,0) 



7(0,1) 



7(M-1,0) 7(M-1,1) 



7(0, TV -1) 
7(1,7V-1) 

7(M-1,7V-1) 



(1.1) 



For the sake of simplicity we will assume that M — N. 

It is also possible to represent the M x N image as a column vector with MN components. 
This can be done by stacking the columns of the M x N matrix into a single vector. This is 
called the vect operation. 

7(0,0) 



i = vect(I) — 



/(M-1,0) 



1(0, N-l) 



I(M - 1,7V - 1) 



(1.2) 
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Figure 1.2: (a) uniform and (b) nonuniform sampling 

We will adopt the following notation. Matrices and vector will be represented with bold letters. 
Matrices will be represented by capital letters. 

Two simple operations we can define on the space of M x N images are the addition of 
two images and the multiplication of an image by a constant (scalar). These operations can be 
easily defined using matrix (left) and vector (right) notations 

W = U + V w = u + v 

(1.3) 

W = aXJ w = au 

The set of M x N images with these two operations is a vector space of dimension MN. 

1.3 Continuous Images 

Continuous images are defined in an infinite set of points. Is it possible to store continuous 
images in a computer? 

At a first glance the answer is negative. A continuous image is defined in an infinite num- 
ber of points and therefore requires an infinite amount of information. This difficulty can be 
overcome if we impose additional restrictions e.g., assuming that the image is an interpolated 
version of a finite set of samples (see figure 1.2a). The problem can be stated as follows. 

Problem: given I(k, I), k, I — 0, M — 1 we wish to compute /(x) where xGl 2 . 

This is known as the interpolation problem [2]. The interpolation is a key step in many 
image processing algorithms which require resampling the discrete image at different locations 
(e.g., rotations, translations with sub pixel resolution). These operations are used in many 
applications e.g., in medical imaging, remote sensing, surveillance, 3D vision, etc. 

The classic solution consists of approximating the continuous function as a sum of interpo- 
lation functions centered at the nodes of a regular grid, multiplied by the image value at those 
points. Therefore, 
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J(x) = J ( k M x - k ) (!- 4 ) 

k£Z 2 

where k G Z 2 and 0 is an interpolation function defined in R 2 . 

The interpolated function I(x) should match the samples values at the nodes of the inter- 
polation grid, /(k). Therefore, the following property applies 

f 1 ifx = 0 

0(x) = { (1.5) 
[ 0 if x G Z 2 \ {0} 

where 0 = (0, 0) is the origin. It is also useful to use separable functions 

0(x) = <f>o(x 1 )<f>o(x 2 ) (1.6) 

where x = (xi, X2). 



Examples 

Some interpolation functions defined in R are 



1 ,/ x 1 - ' l x l if l x l < 1 ,/ x 5in(7rx) , . 
<l>{x) = { 0.5 if \x\ = 0.5 = { <J)(x) = ^ (1.7) 

c.c. 7rx 




These functions verify condition (1.5). This condition is not unique. We can impose ad- 
ditional conditions. For example, given a constant discrete image, /(k) = 1, the interpolated 
function should also be constant /(x) = 1, Vx G M 2 . 

Using (1.4) we obtain 

^0(x-k) = l VxGM 2 (1.8) 

k 

A set of functions (j> verifying this condition is called a partition of the unity. 

The interpolation of a discrete image defined on a regular grid can be performed in a different 
way assuming that /(x) is the sum of base functions ?/>(x) centered at the nodes of a regular 
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grid, multiplied by appropriate constants 



7(x)= ^c(kMx-k) 



(1.9) 



kez 2 



The difference with respect to the previous approach lies on the fact that condition (1.5) may 
not be valid since we do not require that c(k) = ^(k). This approach is known as the generalized 
interpolation problem and allows a greater flexibility in the choice of the basis functions. A 
popular choice for ip are the B-spline polynomials. 

The computation of the coefficients c(k) involves the solution of a linear system of equations 
or filtering the discrete signal /(k) using an appropriate filter. In the case of the B-splines, these 
operations can be done in an exact way using the algorithms described in [1]. 

A more general question can now be asked: can we still reconstruct I(x) when the samples 
are obtained at nonuniform positions? 

This is a more difficult problem but it still has practical interest. One of the motivations 
for this problem is in scope of data compression since we may wish to represent an image by a 
smaller set of chosen samples. 

Most of the information conveyed by an image is located in regions where intensity rapidly 
changes (e.g., close to the edges). To represent an image by a set of points, the sampling 
density should be high at the transitions and sparse at homogeneous regions. The nonuniform 
interpolation algorithms are very useful since they allow to recover the original image from a 
small fraction of its elements. 

An interpolation strategy consists of making /(x) equal to the nearest sample (nearest 
neighbor method) in the image plane x^) 



This is a simple idea. However, the interpolated image is discontinuous and its intensity at 
location x depends on a single sample. 

There are better ways to solve this problem e.g., assuming that the image is band limited non 
uniform reconstruction algorithms can be used. Another alternative consists of approximating 
the image by the superposition of spline basis functions located at the nodes of a regular grid 
[3], or splitting with polygonal meshes. This is an active research topic. 



/(x)=/(x (1) ) 



(1.10) 
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1.4 Color images 



Many images obtained with CCD cameras are color images. This raises a question which we 
should address: 

what is color and how do we represent color images in a computer? 

Color is a human perception. It is the output of the human visual system to electromagnetic 
waves (light) arriving at the retina cells. In the 19th century, Young and Helmholtz found that 
it is possible to synthesize most colors by mixing three primary colors. This is a remarkable 
discovery since it allows us to display color in a computer or in a TV. 

Suppose we receive a light ray at the retina with spectrum 5(A) where A is the wavelength 
and suppose that the three primary colors have known spectra Pj(X),j — 1,2,3. In general, 
5(A) is a complex function and cannot be represented by a linear combination of three spectra 



where cj are the mixing coefficients. 

How do we combine this fact with Young and Helmoltz discovery? Three primary colors 
are not enough to accurately synthesize the input spectrum 5(A). The amazing thing is that, 
despite this difference, the color perception can be the same. 

To understand this we must understand how retina works. Retina has three different types 
of cells sensitive to color, denoted as cones. Two spectra may be different and still produce the 
same response in the three types of cones. 

The output of the cone i G {1,2,3} to the incident spectrum can be modeled as 



where Si(X) is the sensitivity curve of cone i. The sensitive curve has an intuitive meaning: 
Si(Xo) is the output of cone % when it is exposed to a monochromatic wave with spectrum 



Pj(X) i.e., 



3 




(l.ii) 




(1.12) 



tf(A-Ao). 

Using (1.12), we conclude that 5(A) and J2 c jPjW represent the same color iif 1 




i = 1,2,3 



(1.13) 




3 



i = 1,2,3 



(1.14) 



1 iif means if and only if 
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This is a system of three linear equations. If we define, 

Rij = J P j (X)S i (X)d\ (1.15) 

n = J S(X)Si(X)d\ (1.16) 

then 

Rc = r (1.17) 

where R = [Rij],r — [rf],c = [ci\. These are the conditions which must be met in order to 
obtain the same color perception with three primary colors. The mixture is realizable if all the 
coefficients are non negative Q > 0 (additive mixture). This is not true in some cases. When the 
solution of (1.17) leads to negative coefficients the color is not realizable by an additive mixture 
of three primary color sP/. (A) 2 . Can three primary spectra Pfc(A) synthesize all the colors with 
positive coefficients? the answer is no. 

Color can be represented by several systems of coordinates. The most popular is the RGB 
system which is based on three monochromatic colors with wavelengths 700nra (red), 546. lnm 
(green) and 435. Snm (blue). Each color is represented in this system by three coordinates which 
are the mixture coefficients required to synthesize the color. Other systems of coordinates such 
as XYZ, YUV, HSV (hue, saturation, value) and CMY used for printing. 

A color image / is a signal which assigns a 3 color coordinates to each point x G V. Each 
value is therefore a 3-vector. As alternative we can think of a color image as three gray level 
images representing the color coordinates. Figure 1.3 shows an example of a color image and 
the three color components in the RGB system (2nd line) and HSV system. 

Each coordinate conveys useful information. For example, the hue is often used to discrim- 
inate colors and to detect the skin region in an image. 



2 adding more primary colors might help in this case 
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Chapter 2 

Filtering 



2.1 Introduction 

Image filtering aims to eliminate undesired structures from the image and enhance others. Fil- 
tering is used to smooth images, to reduce noise or to enhance intensity transitions or corners. 
We can design filters to enhance vertical or horizontal transitions, squares, circles or to inter- 
polate the image when a few samples are not observed. Figure 2.1 shows a couple of examples. 

In all these cases, the filter maps an input image / into a modified image J. The map can 
be briefly stated as follows 

J = T{I) (2.1) 

There are several types of filters. This sections considers three classes: linear filters, order 
filters and morphological filters. Linear filters are very well studied from a theoretical point of 
view since there is a comprehensive theory to describe their behavior and stability. However, 
non linear filters often achieve better results. 




Figure 2.1: Filtering examples: input images (left) and filtered images (right) 
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Figure 2.2: Linear filtering 



2.2 Linear Filtering 

This section addresses linear space invariant (LSI) filters. A simple example is a filter which 
computes a local average of the input image in a 3 x 3 neighborhood of each pixel (see Fig. 2.2). 
This operation involves multiplying the input image by a group of 9 coefficients which is often 
denoted as mask. The coefficients are moved along the input image as shown in the figure. 

Let us now define two key concepts for this discussion: linearity and space invariance. 

Definition (linearity): A filter T is linear if it obeys the superposition principle 



T(ah + bl 2 ) = aT( Ji) + bT(J 2 ) 



(2.2) 



for all input images h^h and constants a, b. 

This has an intuitive meaning. We can either filter a linear combination of two images or 
filter each them separately and combine the outputs. This property strongly restricts the class 
of transformations allowed. It is interesting to note that a constant map T(I) — C where C is 
a constant image is not linear (unless (7 = 0). 

Linearity is very useful when we decompose the image as a linear combination of basis 
images since we can express the output of the filter as a linear combination of the filtered basis 
images. Specifically, if 

/(k) = 5>0i(k) (2.3) 



then 



^( k ) = ^2 C i^( k ) 



(2.4) 



where ipi — T((j>i). The coefficients remain the same 1 . 

The second restriction comes from space invariance. A space invariant filter is a filter such 
that if the input image is translated by an amount dGl 2 , the output is also translated by d. 



1 For the sake of simplicity we will assume that all the images are discrete images in this chapter. 
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Definition (space invariance): A filter is space invariant if 

T(/(p - d)) = J(p - d) (2.5) 

for all J = T(I) and displacement d G M 2 . 

A filter is linear and space invariant (LSI) if both properties hold. 

We can now state the most important result for LSI filters. If a filter is linear and space 
invariant, its output is given by the convolution sum 

J(k) = ^/(i)H(k-i) (2.6) 

iez 2 

where u(k) is the input image ans h(k) is the response of the filter to a an impulse image S(k) 
(S(k) — 1 if k = (0, 0), S(k) — 0 otherwise). The convolution sum is often denoted by 

J = /*F (2.7) 

The proof is very simple. It is based on the representation of / as a sum of impulses (please 
convince yourself that this is true) 

J(k) = £/(i)«5(k-i) (2.8) 

iez 2 

Applying the filter operator T to both members of this equation 

v(k) = T(u(k)) = T ( £ /(i)J(k - i) ) (2.9) 
Viez 2 / 

v(k) = T (^( k - 1)) linearity (2.10) 

iez 2 

i?(k) = ^ h(k — i) space invariance (2.11) 

iez 2 

as wished. 

This result is remarkable. It states that the filter behavior is completely characterized by 
the response of the filter to an impulse input. To design a LSI filter all we have to do is to 
define the impulse response H(k). 

It is more intuitive to define a set of coefficients (mask) M(k) which is moved along the 
image and multiply the image at each position (see Figure 2.2). The filter mask is the impulse 
response after an horizontal and vertical flip, 

M(m, n) = #(-ra, -n) (2.12) 
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The convolution sum can now be rewritten as follows 

J(k) = £ /(i)M(k - i) (2.13) 

iez 2 

This equation is more intuitive. The filter output is obtained by sliding the mask along the 
input image. Then both images are multiplied and added. Many filtering operations are done 
in this way using (2.13). 



Example 

Consider filter with impulse response H y mask M and an input image / (the origin of each 
image is displayed in bold) 



H 



1 


-1 


M = 


2 


1 


1 = 


1 


2 


1 


2 




-1 


1 




4 


3 



(2.14) 



The filter output is 



J = 




(2.15) 



The properties of LSI filters are well studied. Stability is one of the most important ones. 
A filter is stable if for every bounded input the output is also bounded. It can be easily proved 
that a LSI filter is stable if (and only if) the impulse response is absolutely summable i.e. 

X>(k)l<«> (2.16) 

kez 2 

All the filters with finite impulse response are stable. Only filters with infinite impulse response 
may be unstable. For example, a filter with a constant impulse response h{k) — 1 is unstable. 
One way to define a filters with infinite impulse response is by a recursive equations e.g., 

J(ra, n) — bwl(m, n) + a±o * J(m — 1, n) + aoi * J(m, n — 1) (2.17) 

The equation allows us to filter the input image / with a small number of operations per sample. 
Examples of linear filtering are: 
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Figure 2.3: Gradient vector field 



Image smoothing: a smoothing effect (low pass filtering) can be obtained by computing 
a weighted average of the input pixels. The mask coefficients should be positive and sum 1. A 
typical choice is the Gaussian mask 



|k| z 



M(k) = Ce~^ 



(2.18) 



k = (m,n), ra,n G {— M, ...,M} Partial derivatives: Partial derivatives can only be defined 
in an exact way for continuous images. In the case of discrete images, we can only compute 
approximate values but they are still very useful. Partial derivatives are usually approximated 
by filtering operations 

Hi * / - — & H 2 * I 



dl dl 
— « H ± * I — - 

OXi OX2 

using appropriate masks Mi,M 2 e.g., the Sobel masks 



(2.19) 
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(2.20) 



Image gradient: The image gradient is a 2D vector of partial derivatives 



<ix 



di 

dxi 

dl 

8x2 



(2.21) 



which can be computed using the Sobel masks for example. The norm of the gradient and its 
direction convey useful information about abrupt changes of intensity which can be used for 
edge and corner detection. Figure 2.2 shows the gradient computed for two synthetic images 
using the Sobel masks. It displays the original image and the gradient vector at each point. 
The gradient has a large magnitude close to the transitions and it is orthogonal to the level 
curves of the intensity function. 
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Figure 2.4: Real part of complex exponentials 



2.3 Frequency Analysis 



Frequency analysis is a useful tool for the design of linear filters. 

It is based on the representation of images as a sum of complex exponentials 



(2.22) 



where oo — (001,002) is a vector of frequencies along the coordinate axis. The complex exponen- 
tials are periodic in 001,002 with period 2tt. 



This means that we do not have to consider the whole frequency plane but only an interval 
[— 7r,7r[ 2 since all exponentials are identical to one of these. Figure 2.3 shows the real part 
of complex exponentials for different values of 00 G [— 7r,7r[ 2 : \\oo\\ controls the period of the 
waveform and the direction of 00 defines the direction of the periodic pattern. 

The Fourier transform of a discrete image and the inverse transform are defines as follows: 



e j'(^+27ri)-k _ e J^-k 



Viez 2 



(2.23) 




Fourier transform 



(2.24) 



kez 2 




inverse Fourier transform 



(2.25) 



25 



Figure 2.5: Fourier transform of 8 images: image (left) and magnitude of the Fourier transform 
(right) ) 

where I(cj) is the spectrum of the image /(k). The inverse transform states that the image 
can be synthesized by "adding" complex exponentials with infinitesimal amplitudes I(oj)duj 
(synthesis equation). 

The spectrum is a complex function and periodic with period 2tv as a function of cji and 
It can be characterized by its modulo |/(cj)| and phase arg/(cj). 

Figure 2.3 shows some images (left) and the magnitude of the Fourier spectrum in the in- 
terval [— 7r,7r]. The spectrum of sinusoids with frequency c^o £ [— 7r,7r] 2 are impulses in the 
frequency plane located at c^o- When oo increases the frequency impulses move apart. Further- 
more, if the image undergoes a rotation of amplitude 0 the same happens to the spectrum. The 
second line shows rotated versions of a Gaussian function (Gabor Functions). The spectrum 
belongs to the same family. Is can be noted the effect of scale changes. A compression in space 
brings a dilation in the frequency domain and vice- versa. 

The usefulness of this transform comes from the following properties: 

• The response of a LSI filter to a complex exponential e ja; k is a complex exponential with 
the same frequency and different amplitude 

H(co)e juj (2.26) 

where H(u) is the Fourier transform of the filter impulse response and it is called the 
Frequency Response of the filter and the Fourier transform (analysis equation) tells us 
how to compute the coefficient of each exponential. 

Complex exponentials are very special signals since they are not modified by LSI filter 
T, apart from the amplitude. They are called eigen functions of the linear filters. The 
Fourier transform is therefore an expansion of the input signal as a sum of the filter 
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Figure 2.6: Frequency responses of Gaussian and Sobel filters as a function of ui,U2 
eigenfunctions. 

• The Fourier transform of the filter output is given by the multiplication of two spectra 

J(u) = H(u)I(u) (2.27) 

The Fourier transform of the output image J(u) is the product of the input Fourier 
transform I(uS) by the frequency response of the filter H(k),I(k). 

The second property extends the same idea for multiple exponentials. The input can be 
decomposed as a sum of complex exponentials by the Fourier transform. The amplitude of each 
of exponential is multiplied by the filter frequency response H(cj). 

If |-H"(cj)| is larger than 1 the exponential is amplified, otherwise it is attenuated. We 
can classify filters are low pass, high pass, bandpass, band reject, according to the range of 
frequencies they attenuate. 

Fig 2.3 shows the amplitude of the frequency response for the Gaussian and Sobel filters. 
The frequency response of the Sobel filter is 

Hi(u) = e ~j(-^-^) + 2e~ J ' (_Ct;i) + e - j (- LJl + LJ2) - e -j(-"i-"*) _2e- J ' (_Ct;i) - e -j(-"i+"*) (2.28) 

= (e j(Jl - e - j(Jl )(e j(J2 + 2 + e~ jLJ2 ) = 4jsino;i(l + coscj 2 ) (2.29) 

The frequency response can be used to understand the filter behaviors in the frequency 
domain and can also be used for filter design if we know the ideal frequency response for the 
problem. 
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2.4 Matrix Notation 



We can represent finite signals and images by vectors. Linear filtering operations can therefore 
be expressed as matrix multiplications. 

Let us consider the ID signals and filter first. The convolution of two ID signals v — u*h 
of finite length n, m, can be written as follows 



v = Hu 

where u, v are column vectors and H is a sparse (n + m — 1) x n matrix 



MO) 

ft(l) fc(0) 



(2.30) 



H = 



h(m - 1) 



ft(l) fc(0) 



Ml) 
h(2) 



h(Q) 
h(l) 



(2.31) 



h(m - 1) 

Matrix H is a Toeplitz matrix since the elements of all diagonals are equal: h(i,j) — h(i — j). 
The first column of H is the filter impulse response. 

Let us consider now image filtering. The convolution of two finite support images V — H*U 
with sizes M x M, N x N can also be expressed in the same way 



v = Hu 



(2.32) 



where u, v are again column vectors with all the image elements. However, the structure of H 
is more complex. H it is a block matrix; each block H^- represents the influence of the j-th 
column of the input image on the i-th column of the output. After some manipulations we 
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conclude that 



where 



H = 



H 0 

Hi H 0 
H m 



H*. 



h(0,k) 

h(l,k) h(0,k) 
h(m, k) 



Hi H 0 



Hi H 0 
H 2 Hi 



H 



m — 1 



h(l,k) h(0,k) 



H^ 



ft(0,fe) 
ft(l,jfe) 



H is therefore a block Toeplitz matrix and each block is also a Toeplitz matrix. 
This result can be derived as follows. The column i of the filtered convolution the pixels of 
the i-th column are given by 



(2.33) 



(2.34) 



V (m, i) = ^2^2 H ( m -P,i- 0)U(p, q) 



(2.35) 



p q 



Let us consider the terms which depend on the j-th column of the input q — j and forget 
the others for the moment 



V(ra, i) — H(m — p,i — j)U(p,j) + other 



(2.36) 



This is the expression of a ID convolution of H(m,i — j) * Uj. Therefore the block H^- is 
a Toeplitz matrix whose first column is the impulse response h(m,i — j), ra = 0, • • • , m — 1 as 
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we wished. 
Example 

The convolution of Example 2.14 can be formulated using matrix operations 
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(2.37) 



Although H is block Toeplitz it is not a Toeplitz matrix. Therefore the 2D convolution is 
not a ID convolution. 



2.5 Median Filter 

The median filter is a nonlinear filter. It does not verify the linearity condition. It often performs 
better than linear filters, as we shall see, but it is harder to characterize its performance. 

Median filter is based on the concept of order statistics. Let u\, 1x2, • • • , u n be a set of iid 
random variables, with a common density p(u) and distribution P(u). The order statistics are 
obtained by arranging the n variables in ascending order of magnitude 

< U(2) < • • • < U(n) (2.38) 

is called the k — th order statistic. 
The density of the k — th order statistic is given by 

p k (u) = np(u) I U ~ 1 I Piuf-^l - P(x)) n - k (2.39) 

This can be proved as follows 

Pk(u)du — P{u^ e]u, u + du[} (2.40) 
= nP{u(i) e]u, u + du[ and (k — 1) of the others smaller than u} (2.41) 
= nP{u(i) e]u, u + du[} P{(k — 1) of the others smaller than u} (2.42) 

= np(u)du I n_1 ) P(u) k - 1 (l-P(x)) n - k (2.43) 

Some of the order statistics are well known: is the minimum, is the max- 
imum, and i£( r ), (r = (n — l)/2, n odd) is the median. Consider for example a sequence 
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Figure 2.7: Median filtering: input (up), median filtering (middle) and linear filtering (down) 

{0, 3, 12, 4, 1, 2, 2}. To compute the median, we order the 7 values (0, 1, 2, 2, 3, 4, 12) and pick 
the fourth (u me d — 2). 

The median filter does something similar. It computes the median of the input image in the 
vicinity of each point p 

v(k) = med {u(q) : q G V p } (2.44) 
where V p is the vicinity of p e.g., defined by 

V P = {q:\q-p\<r} (2.45) 

This definition is valid for ID and 2D signals. 

The median filter rejects outliers and preserves discontinuities. Figure 2.5 shows an example. 
The first row displays two input signals: a step signal and an outlier, with and without Gaussian 
noise. The second row shows the output of the median filter of length 5. The median filter 
preserves the transitions and eliminates the outlier as we wish. The Gaussian noise is also 
attenuated. The third row displays the results of a linear filter with triangular impulse response 
of length 7. The linear filter attenuates the Gaussian noise even better but it smooths the 
discontinuity and does not eliminate the outlier. 

The smoothing effect caused by the linear filter is a major drawback in image processing since 
transitions correspond to image edges and they are perceptually very important. Smoothed 
images are perceived as blurred and this may be more annoying than the Gaussian noise itself. 

Another example of median filtering is shown in Fig. 2.1 (left) in which the input image 
is corrupted by "salt and pepper" noise with probability 0.1, i.e., 10% of the pixels were not 
observed and they are randomly replaced by 0 or 255. The performance of the median filter is 
very good in this problem. 
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Figure 2.8: Median filtering: original signal (left), distorted signal with salt and pepper noise 
(center) and filtered signal (right) 



The median value v minimizes the L\ norm of the data 



c = $>fa)-»l 



(2.46) 



qev 



This property is important because it allows several extensions of the median filter. 
Median filtering of vector data 

How can we apply median filtering to a color image? It is not a trivial question because it 
is not possible to order vectors in a meaningful way. We can apply the median filter to each 
component of the vector image but this is not a good idea. Instead, the median filter can be 
defined by the minimization of (2.46) where u(q),v are vectors in this case. This minimization 
cannot be done analytically and have must use numerical optimization methods. 

Weighted median filter 

The median filter assigns an equal weight to all the data points, despite their position inside the 
window V p . It is often useful to assign different weights to the input values according to their 
position. For examples, pixels close to the center of V p are probably more important and should 
have a higher weight the pixels close to the boundary. Weighting can be done by modifying the 
cost function 



qev 

If the weights are integers, the minimization of (2.47) is straightforward. All we have to do is 
to repeat each input data w(q) times and compute the median of the extended sequence. 




(2.47) 
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2.6 Order Filters 



The median filter belongs to a wider class of order filters of the form 

n 

v(k) = Y / a k u (k) (k) (2.48) 

k=l 

where is the kth order statistic of {u(q), q G Vp}, n — H=V P and a k are the filter coefficients. 
A special case are the rank order filters 

v(k)=u (k) (k) (2.49) 

If k — 1 the rank order filter is called the minimum filter and if k — n it is called the maximum 
filter. 

Another order filter is the a-trimmed filter which discards a fraction a of the input samples 
and computes the average value of the others (we assume that a = j/n,j G N) 

«00 = ^4^) Y «(*) (2-50) 

v 7 k=l+a 

The alpha trimmed filter discards the smallest and largest samples. It can therefore deal with 
outliers provided that their percentage is smaller than a. The a trimmed filter performs better 
than the median filter since it is able to reject the outliers and performs a better noise reduction. 

Statistical formulation 

The coefficients can also be obtained from a statistical formulation of the problem [1]. Suppose 
the signal is locally constant and corrupted by additive noise 

u(k) =0 + w(k) (2.51) 

where 9 is a constant and w(k) is a sequence of iid noise variables. We want to find the filter 
with minimum variance (to simplify the notation, we dropped the dependence on p) 

V = E{(6-v) 2 } (2.52) 

assuming that the filter is unbiased 

E{v} = 0 (2.53) 

The bias condition can be written as 

n n n 

E{Y^ a k (6 + w {k) )} =6Y,a k + Yl a kE{w {k) } = 9 (2.54) 

k=l k=l k=l 
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Since the density of is symmetric, the optimal coefficients are symmetric (a& = a n+ i_^) 
and E{w(ty} — —E{w^ n+ i_ k ^}. Therefore the previous equation can be simplified 



53 a* = 1 (2.55) 

k=l 

and the filter output can now be written as follows 

v = 0 + a T w (2.56) 

where a — [a\ . . . a n ] T is the vector of coefficients and w — [w^ . . . w^] T is the vector of noise 
order statistics. The variance is now given by 

V = a T Ra (2.57) 

where R — E{ww T } is the correlation matrix of the noise order statistics. The minimization 
of V with the bias restriction can be obtained by minimizing the Lagrangean function 

L = a T Ra - X(e T a - 1) (2.58) 

where e = [1 . . . 1] T . Computing the derivative with respect to vector a we obtain 

(R + R T )a- Ae = 0 a = ^# -1 e (2.59) 

Since e T a — 1, A = 2/{e T R~ 1 e). The optimal coefficients are 



(«»> 



This method provides an optimal set of coefficients depending on the noise statistics. The only 
difficulty is the computation of the correlation entries Rij — E{w^w^} 



R ii = J W(i)W(j)P(w(i),W(j))dw {i) dw {j) 



(2.61) 



There are analytic expressions for p(w^w^) but the expected value has to be computed by 
numerically integration in the plane. Recent works have try to overcome this difficulty [4]. 
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Chapter 3 

Image Alignment 



3.1 Introduction 

We often have multiple images of the same object and wish to align them in order to compare 
them or to fuse them into a single image. This is what happen in a mosaicing task which aims 
to reconstruct a panoramic image (mosaic) from several partial views. Figure 3.1 shows an 
example in which two images are aligned to create a mosaic. This example is misleading. It 
is not possible to align images of 3D objects by a single transformation unless the objects are 
approximately planar or the camera motion is a pure rotation. This is what happens in this 
example. 

Image alignment is also very important to compare medical images of the same patient 
obtained using different modalities (e.g., CT, MRI) at different instants of time. Alignment is 
also a key operation in surveillance or in stereo vision. 

Our first attempt to formulate the problem is the following: 

Problem 1: alignment of points 

Given a set of points in an image T and a set of points y^ in an image I (see figure 3.2), 
we wish to find a geometric transformation mapping the points into the points y^. 

Image T is called a template and the points x^y^ are denoted as marks or fiducial points. 
Everything becomes more difficult if the matching between the two sets of points x^y^ is 
unknown. 

Alignment can also be done without marks. In this case we may assume that all the points 
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Figure 3.1: Image alignment: pair of images (1st row) and mosaic (2nd row) 




Figure 3.2: Point Alignment 

of the template image T should be associated to points of / with similar color or intensity. This 
is called the color constancy hypothesis. 

Problem 2: dense alignment 

Given two images I,T we wish to find a geometric transformation W(x) which maps points 
from the template image T into points of I, such that J(W(x)) = T(x). 

The alignment of isolated points approach will be discussed in the next sections. Dense 
alignment will be addressed later in this course. 
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3.2 Geometric Transformations 



There are many geometric transformations which can be use to align pairs of images ranging 
from simple translation and rotations which preserve metric information (distances and angles) 
to more general transformations (e.g., projective transformations) which do not preserve metric 
information. 

We will briefly review the most frequent transformations used in image alignment: 



• Translation: 

x + ti 



W(x;0) = 



(3.1) 



where 9 — (ti, t2) T is a translation vector. This transformation preserves distances, angles 
and parallelism of lines. 

• Rigid body: 

W(x;6>) =Rx + t (3.2) 

9 — (R, t) where R G l 2x2 is a rotation matrix and t G I 2 is a translation vector. 
Rotation matrices are unitary matrices (R T R = I) which do not perform axis inver- 
sions (det(R) = 1). The rigid body transformation performs a rotation followed by a 
translation. It preserves distances, angles and parallelism. 

The alignment of two or more sets of points using a rigid body transformation is known 
as Procrustes analysis. 1 . 

• Euclidean similarity: 

W(x;0) = sRx + t (3.3) 

9 — (s, R, t) where s is a scale factor, R G l 2x2 is a rotation matrix tGl 2 a translation 
vector. The similarity transformation performs a rotation followed by a translation. It 
preserves angles and parallelism of lines. 

• Affine transformation: 

W(x;<9)=Ax + t (3.4) 

9 — (A,t) where A G l 2x2 is a square nonsingular 2x2 matrix and t G I 2 is a 
translation vector. This transformation performs rotation, translation, horizontal, vertical 



1 Procrustes was a bandit which offered a stay to tired travelers saying that he had a special bed which would 
fit to the traveler size. He did not explain the way he fitted the bed to the traveler. If the traveler was short 
he/she would be stretched until he/she fit. If it was tall he/she would be amputated. 
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Figure 3.3: Geometric transformation of a square (left) with translation, rigid body, Euclidean 
similarity, affine transformation and projective transformation 

and diagonal scaling. It preserves parallelism of lines but is does not preserve distances 
and angles. 

• Projective transformation (homography) : 

W(x;6>) = 



0ix+0 2 y+03 

6 4 x+6 5 y+6 6 
6 7 x+6 s y+6 9 



(3.5) 



9 — (#i,...,#g). This transformation maps straight lines into straight lines but does 
not preserve distances, angles or parallelism. We note that the vector 9 is defined up to 
a non zero multiplicative factor. The transformation does not change if 9 is multiplied 
by a constant W(x;a#) = W(x;0). This ambiguity consists can be solved by adding a 
additional constraint e.g., ||0|| = 1 where ||.|| stands for the Euclidean norm of a vector 2 . 

All the previous transformation map straight lines into straight lines. There are other 
more flexible transformations which map lines into curves. 



3.3 Aligning sets of points 

Let us assume that we want to align two sets of points x,y G M 2xn . First we must choose the 
geometric transformation W(x, 0). Then we need to estimate the alignment parameters 9. Let 
us address this step. 

The main idea is the following. The parameters should be chosen in order to make the 
alignment errors in the image plane 

e; = y;-W(x i; #) (3.6) 

as small as possible. 

2 the Euclidean norm is denned by ||x|| = V x T x 
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In general we have more constraints than degrees of freedom and it is not possible to make 
all the errors equal to zero. Therefore, a fitness criterion is used to measure the amplitude of 
the alignment errors. The most popular choice is a quadratic criterion (error energy) 

E(0) = , £\\y i -W(x i ;0)\\ 2 (3.7) 

i 

which must be minimized with respect to 9. 

This criterion can be written in a compact way if we define the data matrices x = [xi, 
. . . ? x n ] , y = [yi , . . . , y n ] . e = [ei , . . . , e n ] with dimensions 2 x n. Since 

the error energy is 

E(0) = tr{ee T } (3.9) 

where tr{.} is the trace of a matrix (sum of the elements of the main diagonal). In the sequel, 
we will minimize E(0) for each geometric transformation. We can use (3.7) or (3.9). We will 
use the second. 




Translation 

In this case e = y — x — tl T (1 G l nxl is a vector of ones) and 

E(0) = tr{(y - x - t l T )(y - x - t 1 T ) T } (3.10) 

E{0) = tr{ntt T - 2(y - x)lt T + (y - x)(y - x) T } (3.11) 

To minimize E we compute its derivative with respect to t and make it equal to zero (necessary 
condition). Then 3 

2nt-2(y-x)l = 0 (3.12) 
t = y-x (3.13) 
where x = x l/n,y = y 1/n. The optimal translation is 

t = y-x (3.14) 

where x, y are the mass centers of the two sets of points x^ e y«. The optimal translation shifts 
points x^ in order to align the mass centers of both sets of points. 

3 see the rules of the derivatives with respect to a matrix at the end of this report 
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Affine transform 

In this case 

y = Ax + t T l (3.15) 

The motion parameters A,t are estimated by minimizing (3.9). The estimation of t, A can be 
carried out in two steps. The optimal translation is given by (prove it!) 

t = y- Ax. (3.16) 

The translation vector should align the mass centers. If we now subtract the mean vector from 
the data we obtain a new set of data points with zero mean x, y which are related by the model 

y = Ax (3.17) 

Matrix A is estimated by minimizing the energy 

E = tr{(y - Ax)(y - Ax) T } = tr{yy T - yx T A T - Axy T + Axx T A T } (3.18) 



E = tr{Ryy - 2R p A T + AR^A T } 



(3.19) 



where 



Xii 



E 



%2i 



E^li^H E^li^2i 
YshiXli E^2i^2i 



Ryy = yy T = 



(3.20) 
(3.21) 



E yh E yuhi 
E hiVu E vh 

The derivative of E with respect to A can be easily obtained using the properties of the 
derivative of the trace with respect to a matrix see (10.33,10.37) in Appendix. Therefore, 



and 



— 2Ryx + 2AR, X 



A — R^xl^xx 



0 



(3.22) 



(3.23) 



Equations (3.23,3.16) define the least squares estimate of the affine transform. 
The algorithm can be summarised as follows. 
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Estimation of the Affine Transform 

Given the data matrices x, y, with dimensions 2 x n: 

• center the data at the origin: compute the mean value of the data x, y and subtract them 
from the columns of x e y; The output is x, y. 

• compute the covariance matrices = xx T , ~Ryx = yx T 

• transform matrix A: A = R ya; R~^ 

• translation vector: t = y — Ax 



Rigid Body transformation* 

The estimation of the rigid body transformation is known as procrustes analysis. In this case, 
matrix A is a rotation matrix. The minimization of E is more difficult since it is a constrained 
optimization problem. However, the final result is simple. It is based on the SVD decomposition 
of the cross covariance matrix R ya A 

The rotation matrix R is a 2 x 2 matrix with a single degree of freedom (rotation angle). 
This happens because R has orthogonal columns with unit norm and determinant equal to 1. 
These constraints can be written as follows 5 

RR T = I detA = l (3.24) 

We will use only the first restriction (R T R = /) hoping that the second condition is enforced 
by the data. 

The translation vector is estimated by aligning the mean vectors as before t = y — Rx. 
We can then subtract the mean vector to all the data and estimate the rotation matrix by 
minimizing (3.19) replacing A by R. 

E = tr{Ryy - 2R p R T + RR^R T } (3.25) 

The last term 

tr{RR^R T } = tr{R T RR^} = tr{TL xx } (3.26) 

4 SVD stands for singular value decomposition of a matrix. 

5 a matrix R £ M NxN s called an orthogonal matrix if R T R = RR T = I 
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does not depend on R. Therefore, the minimization of E is equivalent to the maximization of 

tr{K yx K} (3.27) 

with the restriction R T R = I. This optimization problem can be solved by using the La- 
grangean function 

L(R) = tr{R yx H} + A(R T R - /) (3.28) 
where A is a matrix of Lagrange multipliers. A necessary condition for optimality is 

Therefore 

H yx = AR (3.30) 
where R is an orthogonal matrix. Therefore, 

Ry X B% x = ARR T A T = AA T (3.31) 

We will now use the SVD decomposition of TL yx given by 

K yx = UDV T (3.32) 

where U, V are orthogonal matrices and D is a diagonal matrix. Using (3.31) we obtain 

AA T = UD 2 U T (3.33) 

We conclude that A = UDU T . Replacing this in (3.30) we obtain 

UDV T = UDU T R (3.34) 

Therefore, V T = U T R and 

R = UV T (3.35) 

The rotation matrix it the product of two unitary matrices UV T obtained from the SVD 
decomposition of ~R yx 



Procrustes Analysis 

Given the data matrices, x, y, with dimensions 2 x n: 
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• center the data at the origin: compute the average values of x, y and subtract them from 
the columns of x e y; the output matrices are x, y. 

• compute the cross covariance matrix Ryx = y T x and its singular value decomposition 
R~ ~ = UDV T 

• rotation matrix: R = UV T 

• translation vector: t = y — Rx 



Projective Transform 

The projective transform W(x;p) is a nonlinear function of p. Therefore, the minimization of 
the energy (3.7) is a nonlinear optimization problem which does not have an analytical solution. 
If we want to minimize the energy we must resort to numerical optimization methods. 

These difficulties can be overcome by changing the function to be minimized and converting 
it into a quadratic function of the parameters. 

If we define pi = [piP2Pz] T ,P2 = [p^PsPef ,P3 = [PrPsPdf \ e x* = [x* 1 T ], the projective 
transformation can be written as follows 



Vii 



~T 

X- pi 

~T 

X ?P3 



V2i 



~T 

X i P2 

~T 

X ?P3 



These equations allow us to define two algebraic errors (different from (3.6)) 



(3.36) 



en 



x i Pi 



( x fps)2/ 



li 



e 2i = xf p 2 - (xf p 3 )2/ : 



2i 



(3.37) 



and a quadratic cost function 



(3.38) 



This cost function measures the amplitude of the algebraic errors but it has no geometric 
meaning. This error can be written using matrix notation 



e = Mp 



com 



M = 



0 



0 xf 



X„ 



-x[y2i 



0 X^ X^2n 



(3.39) 
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Therefore C — # T M T M#. Since 9 is to be estimated up to a multiplicative factor, we can 
add an additional constraint 9 T 9 — 1. The minimization of C with this constraint can be done 
using the Lagrangean function 

L = 9 T M T M9 - \{9 T 9 - 1) (3.40) 

where A is a Lagrange multiplier. 

A necessary condition to obtain a maximum of L at 9 is obtained by making the derivative 
of L with respect to 9 equal to zero 

2M T M<9 - 2X9 = 0 (3.41) 

M T M<9 = X9 (3.42) 

Therefore, the estimate of 9 must be an eigenvector of M T M. 

There are 9 eigenvectors for M T M. Which one should we choose. It can be easily concluded 
that we should choose the eigenvector associated to the smallest eigenvalue. The cost function 
can be written as follows 

C = 9 T M T M9 = \ 2 9 T 9 = A 2 (3.43) 

The last quality comes from the fact that we are using normalized eigenvectors 9 T 9 — 1. 

The method we have described is very simple. We compute the eigenvalues and eigenvectors 
of M T M and choose the eigenvector associated to the smallest eigenvalue 6 . Figure 3.4 shows 
the alignment of an image obtained with a web camera with a drawing (regular grid). The 
alignment was done using a projective transformation and four marks (the grid corners). We 
have used the minimum number of marks required to estimate this transformation. Some 
matching errors can be noticed in the middle of the grid probably due to the deformation of 
the paper sheet or nonlinearity of the lenses. 

Matching 

The previous methods assume that the marks x,y detected in both images are manually 
matched. How should we proceed when it is time consuming to match the marks? 

There is no efficient algorithm to compute the optimal solution. This is still an open problem. 
To solve this problem in practice we resort to suboptimal methods such as the RANSAC. 

6 the eigendecomposition can be easily done in Matlab using the following command [l,v]=eig(M'M) 
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Figure 3.4: Image alignment: regular grid (left) observed image (center) and aligned images 
using a projective transformation (right) 

The RANSAC method was proposed by Fishier and Bolles [18] performs many (hundreds) 
of tentative matches and chooses the best at the end. 

It starts by randomly choosing the minimal set of points needed to estimate the transform 
and obtain an estimate for the unknown parameters. For example, in the case of the translation 
transform we would chose one point from each image. 

In a second step we transform all the marks of the first image with the estimated transfor- 
mation and check how many of them are close to the marks detected in the second image. This 
is done by computing the distance from the transformed points of the first image T(x^) to the 
nearest points from the second image yj. The number of matches is called the support of the 
model. 

This procedure is repeated many times. The model with highest support is chosen at the 
end. 

Robustness 

The least squares method used in the section is not robust in the presence of outliers i.e., 
wrong observations which are not explained by the model. The estimates produced by the least 
squares methods are severely affected by the outliers since the quadratic cost ||y^ — W(x^, p)|| 2 
assigns a high penalty to these points and these points dominate the optimization procedure. 
A single outlier is enough to jeopardize the estimation procedure. Figure ? shows the least 
squares estimate of a line with one outlier. 

To solve this problem we must use robust estimation methods. The RANSAC method is 
one possible solution for this problem since it is able to ignore a high percentage of outliers 
with a small degradation of quality. Another alternative would be the use of a different cost 
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function 



C = 5>(y;-W( Xi , P )) 



(3.44) 



where p is a penalty function growing slower that the quadratic penalty when the argument 
tends to infinity e.g., p(x) = tan _1 (||x||). 

The criterion main difficulty concerns the minimization of this criterion which can only be 
one using numeric iterative algorithms. A good introduction to robust estimation methods in 
computer vision can be found in [1]. 

3.4 Alignment without marks* 

This section is based in [2]. We wish to align two or more images without using marks. In 
this case geometric criteria are no longer valid since we do not know any marks. A different 
criterion is used instead based on a color constancy hypothesis. We assume that homologue 
points have similar intensity (color). 

This can be done using a quadratic criterion 



The geometric transformation between a pair of images W(x;0) is estimated minimizing E. 

The assumptions we are making are strong in some applications. Intensity (color) constancy 
is not true in practice if the camera or the object undergoes a large movement. Second the 
motion model is usually not valid for all the image points. In this case, the method can still by 
applied in small regions e.g., 5x5 blocks. 

The minimization of E is a non linear problem. It is therefore necessary to adopt recursive 
numeric algorithms to obtain the motion model. In the sequel we describe the well known 
Lucas-Kanade method. 

Lucas-Kanade method assumes that the current estimate of 0 is slightly modified 



If AO is small we can approximate J(W(x, 0 + AO) by the first terms of the Taylor series 




(3.45) 



X 



9 <- 0 + A9 



(3.46) 



J(W(: 



x; 9 + A0) = J(W(x; 6)) + V/(W(x; 6)) 



T <9W(x;6>) 



A6» 



(3.47) 



Replacing this in the energy expression we obtain 




X 



(3.48) 



47 



Lucas-Kanade Algorithm 

iterate until 9 converges 

• transform / through W(x;0) to obtain 7(W(x, 0)) 

• compute the error image T(x) — J(W(x; 0) 

• transform gradient V through W(x;0) to obtain V/(W(x;0)) 

• compute the images V/(W(x; 9)) T 

• compute matrix R 

• compute vector r 

• compute A0 - R _1 r and update 9 <- 0 + AO 



T r 

T(x) - J(W(x; 0)) - V/(W) T ^ A0 



= 0 (3.49) 



Figure 3.5: Lucas-Kanade algorithm 

The optimization of E with respect to A0 is easy now. Computing the derivative and making 
it equal to zero leads to 

E(™<w,^) 

This is a system of linear equations 

RAi9 = r (3.50) 

where 

B = £^ T WCW)W(W)^ (3-51) 

X 

r = E^ V/(W)(T(x)-/(W(x;0)) (3.52) 

X 

Lukas Kanade algorithm is summarized in figure 3.5. Figure 3.6 shows several steps involved 
in the update of 9 through Lucas-Kanade method. This figure was extracted from [2]. 

Some comments are required. We have assumed until now that /, VI are continuous images. 
The steps involving the warping of the image must be performed using interpolation algorithms 
with sub pixel accuracy. The recursion is computational demanding since the warping opera- 
tions associated to the transformation W(x, 0) must be performed at each iteration. There are 
alternative approaches which are computationally more efficient (see [2]). 
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Figure 3.6: Lucas Kanade algorithm: operations involved in each iteration. Extracted from [2]. 



49 



Let us see a simple example. 



Example - translation 

Suppose the motion model is a simple translation W(x, 0) — x + t, (p = t). The derivative 
of W with respect to 9 is the identity matrix. 




(3.53) 



In each iteration the translation is estimated solving a system of 2 equations Rt = r where 

R = W(W)W(W) T r = ^ VJ(W)(T(x) - J(W(x; p)) (3.54) 

X X 



3.5 Discussion 

The alignment methods based on dense sets of points do not require point matching which is 
usually error prone. There are however other difficulties: a more difficult optimization and lack 
of robustness with respect to outliers and initialization. 

The cost functional is a non convex function and has several local minima. The estimates 
obtained by numeric optimization algorithms get stuck in local minima. Therefore, the final 
estimate depends on the initialization. How can we solve this difficulty? the use of multiple 
scales has been an effective tool to deal with this problem [3] . A first estimate of the parameters 
is obtained using a coarse representation of the image. This estimate is refined using higher 
resolution images. 

The estimation performed by the Lucas-Kanade method is sensitive to outliers. There are 
image regions where matching is not possible due to lack of information or because the model 
is no longer valid and must be replaced by another model. These regions significantly degrade 
the estimates since we are using a quadratic criterion. Robust estimation methods have been 
proposed (e.g., see [4]). 
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Chapter 4 



Image Representation and 



Models 



4.1 Multi-scale representations 

An image contains objects of different sizes e.g., large buildings and smaller people. Small 
objects are best characterized if the camera sticks close to them while large objects should be 
observed at mutch larger distances. This means that we should characterize different objects 
at different scales, the scale being related to the distance from the camera to the object. 

What is a multi-scale representation? a multi scale representations of an image / is a 
description of I which allows the analysis of the image at different scales, ranging from coarse 
scales to fine ones. An example of such a representation is the scale-space concept and the 
Gaussian pyramid considered in the next sections. 

4.2 Scale-Space and Gaussian Pyramid 

Let I(x) be continuous image and let us define a distorted image X(x; s) obtained by smoothing 
the image / with a Gaussian kernel. 



I(x;s) =/(x)*G a (x) 



(4.1) 



G s (x) = 



1 



(4.2) 



e 



2s 
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Figure 4.1: Scale space images for s — 0, 2, 4, 8, 16 



The convolution destroys the details smaller than s. The parameter s is therefore a scale 
parameter and the function X(x; s) is called the scale-space representation of /. 

Figure 4.1 shows some slices of X(x; s) for different values of s showing the way small details 
are being lost. 

When we work at coarse scales (large s) we need less information to represent the image. 
We can therefore devise ways of reducing the amount of information at coarse scales. This can 
be done by making a change of variable. 



If the size of X(x; s) is T x T the size of J(~x; s) for a given s is j x j. This is called a pyramid. 
The use of multiple scales improves the performance of the image analysis algorithms which 
can be tailored for a specific scale and reduces the computation effort since coarser scales can 
be processed much faster. 

In practice we use discrete images and discrete scales values e.g., s — 2 J . The Gaussian 
pyramid can be defined by low pass filtering the image and down sampling by a factor of 2 to 
reduce the size of the image to one half at each iteration. Let / be a N x N discrete image. We 
assume that N — 2 M for the sake of simplicity. We initialize the level 0 of the pyramid as 



and recursively compute the other levels by smoothing the previous level and resampling the 



J(x; s) = X(sx; s) 



(4.3) 



/o(k)=/(k) 



(4.4) 
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Figure 4.2: Gaussian Pyramid, level 0, 1, 2, 3, 4 



smoothed image using with a 2 : 1 factor 

Ji (k) = U (k) * G a (k) (k) = Ji (2k) i = 1 , . . . , M (4.5) 

The sequence of images Iq,I u ...,I m of sizes 2 M x 2 M , 2 M ~ 1 x 2 M_1 , . ..,1x1 define the 
Gaussian pyramid for the discrete image /. Figure 4.2 shows 4 levels of the Gaussian pyramid. 
The original image is shown on the left. 

The Gaussian pyramid has one drawback. It increases the amount of memory needed to 
store the image. This leads to the following question: can we build a multi-scale representa- 
tion without increasing the number of coefficients? This question will be answered using Haar 
functions and wavelet theory. 



4.3 Haar Representation 

Let us consider a ID approximation problem first. Suppose we want to approximate a ID 
function / : R — > R by a piecewise constant function f 0 such that fo is constant between any 
two consecutive integers [k, k + 1[. We can write this approximation as 

/o(a;)= Yl c ^ x ) ( 4 - 6 ) 

k= — oo 

where (j>(x) is the box function (see Figure 4.3 left) 

f 1 0 < x < 1 
<K*) = { ~ (4.7) 

[ 0 otherwise 

and Ck is the average value of / in the fc-th interval [fc, k + 1[. 

We can improve the accuracy of the piecewise approximation by considering smaller intervals 
e.g., by splitting each interval into 2 J subintervals. Every time we split the intervals we are 
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Figure 4.3: Haar basis functions: a) Haar scaling function; b) Haar wavelet 



changing the scale. If we repeat this procedure we obtain a multi-scale approximation of / with 
an accuracy as high as desired. The index j is the scale level or resolution. The approximation 
obtained at level j is 1 

/;(*)= E CjkWx-k) (4.8) 



The basis functions are box functions compressed by a factor of 2 J i.e., with a shorter duration 
and shifted. 

Figure 4.4 shows the multi-scale approximation of a speech signal by piecewise constant 
functions for j — 3, 5, 7 and figure 4.5 shows the evolution of the norm of the approximation 
error when the number of levels j increases 



D(j) = 



e 5 = f ~ fj 



(4.9) 



The error tends to zero as expected. 

The representation (??) is very useful but it has a drawback. It is not recursive. Every time 
we want to increase the resolution level j we have to define a new set of basis functions and 
to compute all the coefficients again. Can we avoid this? The answer is yes and it is very 
simple in this case. 

Let us consider the simplest case of a constant approximation fo(x) in the interval [0, 1[. If 
we want to split this interval in two we only have to consider a new basis function 



ip(x) — < 



1 0 < x < \ 
-1 \ < x < 1 
0 otherwise 



(4.10) 



and multiply it by an appropriate coefficient: fi(x) = fo(x) + dtp(x). This idea applies if we 
want to compute fj + i from fj in K. We just have do add a new set of base function 



Cj = {<p(2 j x - jfe)} 



(4.11) 



Hhe functions are usually multiplied by a gain 2-? to keep their energy equal to 1. We have omitted the 
normalization factor for the sake of simplicity. 
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Figure 4.4: Hierarchical approximation of a speech signal using piecewise linear functions (j — 
3,5,7) 



Therefore, 



+00 



fj+i( x ) - fj( x ) + dj( x ) dj - ^2 d jk xl)(2 3 x - k) 



(4.12) 



k= — 00 



The signal fj is called the approximation of level j and dj is called the detail. The functions 
ip(2 l x — k) are called the Haar functions and ip(x) is called the Haar wavelet. It can be shown 
that they are orthogonal and they span the L 2 (M) space, the space of functions with finite 
energy. Therefore the approximation function fj — > f when j tends to infinity. 
Combining these equations we obtain 



j-i +00 

i=0 k= — oo 



(4.13) 



The function f 0 can also be expressed in terms of wavelets if the scale variable % starts at —00. 
We obtain two alternative ways to express fj using scaled versions of (j> and ip 

+00 



fj( x ) = ^2 c jk (j)(2 J x - k) 



k= — oo 



j-l +00 



(4.14) 



(4.15) 



%— — 00 k— — 00 



This decomposition can be interpreted in terms of subspaces. The function fj is a linear 
combination of linearly independent functions (j){2^ x — k). This means that fj belongs to a 
subspace Vj span by the base 

Bj = {(j)(2 j x -k),ke Z} (4.16) 
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Figure 4.5: Norm of the approximation error as a function of the number of levels 
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Figure 4.6: Slots associated to hierarchical basis functions (3 levels) 



The subspaces Vj are nested i.e., 

V-i C V 0 CV ± CV 2 ... (4.17) 

The subspaces of lower resolution are contained in the subspaces of higher resolution. Since 

fj+i = fj + dj (4.18) 

where the detail dj is a signal which belongs to the subspace Wj — span{i/j(2^ x — k)}. Each 
element of V}+i is a sum of an element of Vj with an element of Wj. Therefore, 

V j+1 = Vj © Wj (4.19) 

where 0 is the direct sum of subspaces. 

What is the relationship between the coefficients of both decompositions (4.14), 
(4.15)? The two sets of coefficients are closely related. If we know the coefficients Cj we can 
compute the coefficients of the approximation and detail of lower levels as follows 

C jk = « (Cj+l,2fc + Cj+ lt 2k + l) djk = x(Cj+l,2fc - Cj+ lt 2k+l) (4.20) 
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For example given the coefficients c^k — (46553324) we can easily compute the approximation 
and detail coefficients for each level by computing sums and differences 

Level 3 46553324 

S \ aproximagao 
Level 2 5533-1 00-1 

, ,h Mdetalhe 
Level 1 5 3 0 0 LJ 

Level 0 4 1 

4 1 0 0- 1 00-1 

coo = 4,d 0 o = Mi* = (00),d 2k = (-100-1). 

These ideas con be extended to a broader class of hierarchical representations as we will 
now discuss. 



4.4 Wavelet Transform 

4.4.1 Multiresolution Approximation 

The Wavelet transform aims to approximate signals at different resolution levels, using scaled 
and shifted versions of a unique basis function called mother wavelet. Figure 4.4.1 shows an 
image (left) and its decomposition using 2 levels of the wavelet transform (right). The sub 
images represent the lower resolution approximation (up/left) and the six detail images i.e., 
information needed to reconstruct the original image. The process can be repeated. The low 
resolution image can decomposed again and replaced by four sub images. This allows to create 
a multi resolution representation which keeps the amount of coefficients constant. The wavelet 
theory is closely related to several topics in mathematics and signal processing: multi-resolution 
approximation, functional analysis and multi-rate filter banks. Excellent accounts on wavelets 
can be found in [1, 2, 3, 4]. 

Let us suppose we want to approximate a function f(x) G L 2 (R) where L 2 (M) is the space of 
functions with finite energy. We will consider a set of nested subspaces of increasing resolution 
{Vi,ieZ} 

...V-i C V 0 CV 1 CV 2 ... (4.21) 

such that if a function f{x) G Vj then f(2x) G Vj+i- This means that if a function f{x) belongs 
to Vj the function obtained by compressing it by a factor of 2 belongs to the next subspace. 
Furthermore we will assume that the limiting set 

lim Vj (4.22) 
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Figure 4.7: Wavelet transform: approximation (upper left block) and 6 details 

is dense in L 2 (M) i.e., any function / G L 2 (R) can be approximated as well we wish by a 
function fj G Vj provided that we choose j high enough. 

Since the subspaces Vj are nested we can write each of them as a direct sum of subspaces 



where V}_i is the lower resolution space and Wj-i is a subspace of detail. Therefore each 
function fj G Vj can be decomposed in a unique way as 



where fj-\ G V}_i and dj-i G Wj-\. These functions are called the approximation and detail 
at resolution j — 1. 

How can we compute the approximation and the detail at several resolutions? 

The answer is easy if we know an orthogonal bases (j>jk(x),ipjk(x) for the subspaces Vj,Wj. 
Since the subspaces Vj at different scales are related by dilations by power of 2 it would be nice 
if the basis functions at different scales were also related by scaling and shifting operations. 
This problem is answered by the next theorem. 
Theorem [Mallat]: 

Given a multi resolution approximation {Vj,j G Z} there exists two unique functions 



Vj = Vj- ± e Wj-! 



(4.23) 



fj - fj-i + 



(4.24) 



(j)(x) G L 2 (R) 



\l){x) G L 2 (R) 



(4.25) 
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called scaling function and mother wavelet such that \2S(\)(2?x — k)} is an orthonormal 
base for Vj and \2h^(2?x — k)} is an orthonormal base for Wj. 
For the sake of simplicity we will use the following notation 

cl> jk (x) = 2^(2 j x-k) (4.26) 

il> jk (x) = 2ty(2ix-k) (4.27) 

for the basis functions of Vj,Wj and the functions ipij(x) are called wavelets. 

If we have orthonormal bases for the subspaces, the projections of a function / on Vj and 
Wj are given by 

fj( x ) = ^2 c j,k^j,k(x) c jik =< f,<f> jlk > (4.28) 

k 

d j( x ) = ^2 d j,k^j,k(x) d jik =< f,il>j,k > (4.29) 

k 

4.4.2 Fast Wavelet Transform 

Equations (4.28,4.29) are not practical since they involve the computation of an integral for each 
coefficient at each scale. There are faster ways to compute the coefficients using the so-called 
fast wavelet transform. 

We have not used yet the fact that the base functions at different scales are related by 
a contraction. This leads to important properties which are at the heart of the fast wavelet 
transform algorithm. The scaling function and the mother wavelet belong to Vb and Wo 7 
respectively. Therefore they are also members of the higher resolution subspace V\ and can be 
expressed in terms of its basis functions 



H x ) — ^2k ^o(fc)0i,ife(aO (self-similarity) 



(4.30) 



where 

h 0 (k) =< (t>(x),(t> lik (x) >= V2 J cj)(x)(j)(2x - k)dx (4.31) 

hi(k) =< ip(x),(t>i j k(x) >= V2 J ^{x)4>(2x - k)dx (4.32) 

Using this result, it is easy to show that the approximation and detail coefficients at con- 
secutive levels are related by (see Appendix) 

c J-hk = X^°( n 2k ) c jn d j-i,k = y^/ii(n - 2k)c jn (4.33) 
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Figure 4.8: Fast wavelet transform: multi rate filter bank with perfect reconstruction 

cjn = y^oQi 2k)cj-i, k + y^fti(n - 2k)dj- 1:k (4.34) 

The first equation computes the pyramid from the high resolution coefficients i.e., it computes 
successive coarse approximations and details of the signal from the high resolution coefficients. 
The second equation does the opposite. It computes the high resolution approximation from 
the low resolution coefficients and details. 

We can easily express these equations in terms of convolution, up sampling and down sam- 
pling operations defined as follows 



c(n/2) n even 
0 n odd 



[t]c(n) : 

The previous equations can be written as 



[i]c(n) = c(2n) 



(4.35) 



cj-i,k = [i](h 0 (-k) * cjk) dj-i,k = [i](hi(-k) * cjk) (4.36) 

cjk = M*0 * [t]cj-i,fc + fti(n) * [t]dj-i,k (4.37) 

These equations define a multi rate filter bank which implements both steps (see Figure 
4.8). This filter bank comprises two parts: the analysis and the synthesis which correspond 
to the computation of the lower resolution coefficients from high resolution coefficients and 
the opposite. The wavelet transform maps a set of high resolution coefficients into a set of 
coefficients corresponding to different resolution level. 

This established a close link between wavelet theory and multi rate filter bank theory. The 
wavelet decomposition is associated with a perfect reconstruction multi rate filter bank. 



4.4.3 Representation of discrete images 

How do we apply the wavelet decomposition to discrete signals? 



62 



The wavelet decomposition provides a hierarchical representation of continuous signals. 
However we can use it to represent discrete signals as well. This can be done by initializ- 
ing the coefficients of the j level with samples of the discrete signal, c J?n = f{n), and then 
apply the fast wavelet transform to obtain the approximations and details at lower resolution 
levels. 

How can we use this transform with images? This question can be solved in an easy 
way if we use separable basis functions. This amounts to computing the wavelet transform to 
each column of the image and the apply it again to each row of the transformed image. 

4.5 Probabilistic Models 

Probabilistic models are used in many image processing operations. Suppose we wish to separate 
objects of interest from a background image. To solve this problem we need a statistical model 
of the background image, in order to detect regions which do not fit the background model. 
This techniques is used in background subtraction methods. Other examples are provided by 
image denoising and restoration problems in which we want to improve the quality of images 
corrupted by noise. Suppose that the observed image is corrupted by additive noise 

J = I + W (4.38) 

in which / is the original image, J the observed image and W is the noise. We need probabilistic 
models for the original image / and noise W in order to separate them. Without models we 
can not separate / from W since we only know their sum. 

The basic question is the following: given a pair of images we want to know which one has 
higher probability, in a certain context. Figure 4.9 shows two binary images which represent 
moving objects in a scene. Which one is more probable? The answer is very easy in this case. 
The left image is the most probable. We can recognize the objects from their silluette. Prior 
information about the objects shape from several viewpoints is an important cue. A simpler 
criterion is that we have chosen the smoother image i.e., the image with less transitions. Can 
we include this kind of evaluation in a probability measure? 

To characterize a random image / = (7i,...,/ n ) we need to know its joint probability 
distribution 

p(I)=p(I 1 ,...,I n ) (4.39) 

which depends on all its (thousands) variables. It is not possible to estimate this distribution 
from experimental data except if we make very strong hypothesis. This can be done if we are 
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Figure 4.9: Are all the images equally probable? 



working with a specific type of images (e.g., normalized faces). However, in the general case it 
is not possible to obtain accurate models in this way. Therefore, the probabilistic models used 
in image processing are poor models in the sense that they only try to represent a small subset 
of the properties of real images. Still, they play an important role in many applications. 

4.6 First and second order histograms 

The histogram is a simple way to characterize an image. To build a histogram we split the 
intensity range into a number of bins and count the number of pixels which fit in each bin. In 
this way we are estimating the probability distribution of isolated pixels p(h) assuming that 
this distribution is valid for the whole image. Figure 4.10 shows the histograms associated to 
two images. 

The histogram is a very crude representation. It neglects all spatial information and it also 
assumes stationarity. If the pixels were independent variables, the joint distribution could be 
computed multiplying these factors 

k 

p(I) = l[pi(Ii) (4.40) 

i=i 

However this is not true. Images are highly correlated. A second step consists of studying pairs 
of pixels with a given spatial relationship e.g., neighboring pixels along horizontal lines. We 
can now define a scatter plot ploting all pairs of image intensities If, Ij such that i,j obey the 
geometric restriction. Figure 4.11 shows the scatter plots for the two previous examples. 

We can characterize pairs of pixels using a 2D histogram which estimates the joint distribu- 
tion of neighboring pixels p(Ii, Ij). The relative position of pixels i,j is defined by the user e.g., 
pixels separated by a distance d along a given direction. The histogram is computed using all 
pairs of pixels with the same relative position belonging to the image or to a specific region of 
the image. We are therefore assuming that the image is a stationary process i.e., its statistical 
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Figure 4.10: images (up) and intensity histograms (down) 




Figure 4.11: Scatter plots and 2d histograms 

properties remain the same everywhere. Of course this is not always true. It is true if we are 
dealing with an aerial image of a forest but not true if we apply the method to the image of a 
face for example. 

The 2D histogram of pairs of pixels is a way to characterize textures in a given region and 
has been used for texture recognition. The 2D histogram is often known as co-occurrence matrix 
in the context of texture analysis. 

Unfortunately the application of 2D histograms is very restrictive and the idea can not be 
extended to characterize the whole image involving thousands of variables. How should we 
proceed then? 
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Figure 4.12: Scatter plots of detail coefficients for the boat (up) and noise (down) images: 
horizontal detail (left), vertical detail (center) diagonal detail (right). The plots were built 
using pairs of consecutive coefficients along the horizontal direction. 

4.7 Wavelet models 

Neighboring pixels usually have a strong dependence. Modeling this dependence is a difficult 
step in the design of a probabilistic model for the image. However, it has been found that the 
wavelet transform decorrelates the image since the detail coefficients are less correlated. 

Figure 4.12 shows the scatter plot for the three details of the first level decomposition of 
images boat and horse?. We have considered pairs of neighbor coefficients along the horizontal 
direction for the three types of detail images. We can observe a strong concentration of points 
in the vicinity of the origin (many coefficients have a small amplitude). There is no apparent 
dependence between the coefficients and it is reasonable to assume that the coefficients of each 
detail image are approximately independent. The same can be said as a first approximation 
if we compare the coefficients of different detail images within the same scale or at different 
scales. If we stack all the detail coefficients in a single vector d — (di, . . . , d m ) we can assume 
that 

n 

P(d)=j[p(di) (4.41) 
Since there is a 1-1 linear map between the detail coefficients and the image I 

d = WI (4.42) 

with det(W) = 1 then 

p(I) = p(d) (4.43) 
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Figure 4.13: Image graphs: 4 and 8 neighborhood systems and corresponding cliques 

4.8 Markov Random Fields 

An image is a set of random variables organized in rows and columns (see Fig ?). There is 
dependence among these variables: if we know the value of some pixels we can predict the value 
of others. In fact there is a short range dependence and long range dependence mechanisms. 

It is perhaps reasonable to assume that the influence of the whole image on a given pixel Ii 
can be discarded if we know its neighbors. This can be written as 

p(Ii/I - {/<}) = pih/Ni) (4.44) 

where I — {Ii} is the set of all image pixels except % and Ni are the neighboring pixels of pixel 
i. This is called the Markov property. 

Let Q — (5, N) be a graph consisting of a set of sites (nodes) S and a set of links N — 
{(hj)jhj £ 3} The sites represent the pixel location and the links represent neighborhood 
connections. Figure 4.11 shows two graphs corresponding to two neighborhood systems: 4- 
neighborhood and 8-neighborhood. A group of nodes fully connected i.e., such that any pair 
of nodes is connected is called a clique. Figure 4.11 also shows the different types of cliques 
associated to the graph. 

Let / = (If, i G S) be a set of random variables associated to the nodes of the graph. In our 
case / is an image and the nodes are the image pixels. 

Definition [Markov Random Field] / is a Markov Random Field (MRF) iif (if and only if) 

i) P(I) > 0, VI 

ii) p(Ii/I - {Ii}) = p(Ii/Ni), Vi G S (Markov property) 

The first is a technical assumption. The second is the Markov property which means that 
the vicinity of a pixel isolates the pixel from the rest of the image. To predict the value of a 
pixel we just have to know its neighbors. 

This seems to be a reasonable assumption which takes pixel dependence into account as we 
wished but builds long range dependences using local dependences. 
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4.8.1 Gibbs Random Fields 

Let us consider a related concept. A Gibbs random field (GRF) is a set of random variables 
/ = G S) associated to a graph Q — (5, N) such that p(I) is a Gibbs distribution 

p(J) = \e~ E ^ (4.45) 

where E(I) is an energy function defined as follows 

E = J2 Vc(Ie) (4.46) 

where C is the set of cliques of the graph, V c is a function of the clique variables I c called poten- 
tial function and Z is a normalization constant called the partition function. This expression 
states that the energy is a sum of potential functions V C (I C ) associated to the cliques of the 
graph. 

The clique potentials V c are specified by the user or estimated from the data. Estimation is 
usually a difficult step in this type of models. The partition function is given by 

z = Y, e ~ E{I) ( 4 - 47 ) 

It involves the sum of the exponential function for all possible images. There is no analytical 
expression for Z in most problems and it is impossible to evaluate it in a computer due to the 
huge number of terms involved in the sum. In most of the cases we do not know Z. 

The probability distribution can be factorized as a product of terms each of them associated 
to a clique of the graph 

p(T) = J] *c(!c) (4-48) 

where 

$ c (/ c )oce- Vc ^) (4.49) 



Example - Ising model 

The Ising model was proposed in the context of statistical physics for the study of ferro- 

magnetism 2 . It assumes that the image / is a binary random field with a Gibbs distribution 

2 Ernst Ising (1998-2000) was a German physicist who studied a model for ferromagnetism in is Ph.D. thesis. 
Ising was persecuted by the nazis during World War II and moved to the US after the end of the war. He taught 
physics and mathematics at Bladley university until his retirement. 



68 




Figure 4.14: Realizations of the Ising MRF (/? = 0.8) 



with energy 

E(I)=/3 £ Kh-Ij) (4.50) 

(ij)eN 

where 5 is the impulse and N is the set of all 4-neighbor pixels. The clique potentials are 
defined as 

V{i,j) = p8(Ji-Ij) ( 4 - 51 ) 

for all neighboring pairs of pixels. In the Ising model, the energy counts the number of horizontal 
and vertical transitions. A constant image has energy E — 0. The parameter (3 is defined as 
the inverse of the temperature and controls the average number of transitions in an image. 
Let us consider an example 

0 0 0 0 

0 110 

0 0 0 1 

0 0 11 

The energy is E — 10/3. 

Figure 4.14 shows images produced by this model for three values of (3 (pixels marked with 
1 are displayed in black). The size the regions depends on /?. This model is very simple but it 
is useful in several applications e.g., to represent moving regions on a background image or two 
types of land cover in aerial images. 
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4.8.2 Equivalence of MRF and GRF 

Hammersley- Clifford proved that Markov random fields and Gibbs fields are equivalent i.e., a 
random image is a MRF on a graph Q iif it has a Gibbs distribution 

p(J) = \e~ E W (4.52) 

where 

£7=X)V C (J C ) (4.53) 
and C is the set of all cliques. This is called the Hammersley- Clifford theorem. 

4.8.3 Simulation 

It is often useful to simulate the MRF and generate realizations of the field variables using the 
Gibbs distribution. Two methods which are often used to simulate MRFs are the metropolis 
algorithm and the Gibbs sampler. 

The metropolis algorithm produces a sequence of images as follows. Given an image / we 
introduce a random change 

(4.54) 

and compute the energy difference 

AE = E(I') - E(I) (4.55) 

If AE is negative the new image V is accepted, otherwise V is randomly accepted with proba- 
bility P — e AE between 0 and 1. If P is close to 1 we accept I' with a high probability. If P is 
close to zero we reject V most of the time and keep the previous image /. The random change 
can be done in many ways. We can make a local change (e.g., modify a single pixel) or a more 
global one (modify a group of pixels). The best strategy depends on the problem. However, it 
can be shown that under broad conditions the sequence of variables /, generated in this way, 
has an asymptotic Gibbs distribution (4.52). 

To generate a realization of the MRF with distribution (4.52) we should perform a large 
number of iterations of the Metropolis algorithm and pick the last image. If we want to obtain 
several realizations of the random field we should not pick the images at consecutive iterations 
since they will be dependent. To obtain independent realizations we must perform several 
iterations between two consecutive realizations. 
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Metropolis Algorithm 


Gibbs Sampler 


initialization: initialize / randomly 


initialization: initialize / randomly 


iterate N times 


iterate N times 


random change: I ^ I 1 


iterate for all pixels 


energy variation: AE — E(I') — E(I) 


randomly change pixel i: p(Ii/Ii) 


decision: accept V with prob. 


end 


P = min(l,e AE ) 


end 


end 


MRF realization: / 


MRF realization: / 





Table 4.1: Metropolis algorithm (left) and Gibbs sampler (right) 



An alternative approach is the Gibbs sampler. The Gibbs sampler starts from an random 
initialization of the image and modifies one pixel at a time scanning the image by rows or 
columns. Each pixel is randomly modified keeping all the other constant. This is done by 
sampling the conditional distribution p(Ii/Ii) where Ii — I — {h}. The conditional distribution 
is also a Gibbs distribution 

pih/ii) = Y^ Ci Vc(I) ( 4 - 56 ) 

where C\ is the set of cliques which contain the pixel i. The only variable in this expression 
which is allowed to vary is U since all the other are considered as constant and remain invariant. 
The partition function Z 1 can now be easily computed 

Z' = ^2 e E ^ Vc(/) (4.57) 
h 

the sum involves a single variable. 

The Gibbs sampler is summarised in Table ??. The convergence of the Gibbs sampler to 
the asymptotic distribution is usually faster that the convergence of the Metropolis algorithm 
but both methods are used in practice. 

4.8.4 Optimization 

Many problems in computer vision and image processing can be modeled by Markov random 
fields with Gibbs distribution 

p(I\J) = ±e~ E ^ (4.58) 
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where J is the available information (observations). We often want to know what is the most 
probable image /. This can be done by minimizing the energy 



/ = argmin E(I; J) 



(4.59) 



In many cases a 4-neighborhood system is adopted and the energy function has two types of 
terms: terms associated to isolated pixels and pairs of pixels 



The minimization of this energy is still a very difficult problem since the energy depends on 
a huge number of variables (typically hundreds of thousands) and it is a non convex function. 
Most optimization algorithms perform poorly in this type of problems and usually get stuck in 
local minima. 

The optimization step has been a major bottleneck which prevents a widespread use of MRF 
in computer vision. Several techniques have been proposed to deal with this problem. We will 
briefly mention a few of them: 

• Iterated Conditional Model (ICM) 

ICM is a very simple method. It basically optimizes one variable at a time keeping all 
the other constant. For example, if / is a binary image we visit each pixel at a time 
and choose the label (0 or 1) which minimizes the energy, freezing all the others. After 
visiting all the pixels we repeat the whole procedure starting from the image obtained 
in the previous iteration. The ICM method is simple but it usually gets stuck in local 
minima and may provide poor results. 

• Stochastic methods 

The idea is also simple. Let us define an auxiliary MRF with Gibbs distribution 



If we increase (3 the distribution tends to zero except in the vicinity of the configurations 
which correspond to the global maximum (see Figure). 




(4.60) 



v (j) - I e -/3£(^) 

Faux y 1 ) — ^ c 



(4.61) 
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Evolution of the Gibbs distribution when f3 increases. 

To obtain the global maximum the MRF is simulated using the Metropolis algorithm or 
the Gibbs sampler. The parameter f3 is slowly increased during the simulation. It can 
be shown that if we increase /? slow enough the sequence of MRF generated in this way 
converges to the global maximum. These methods perform well but they are very slow. 

• Optimization on graphs 

It has been recently proposed a class of optimization algorithms based on graphs. These 
methods perform better than the previous ones but their theory is still a research issue 
and it is outside the scope of this text. The interested reader may find an introduction 
to this topic in [?]. 

4.8.5 Parameter Estimation 

The energy function often depends on parameters chosen by the user. For example, the Ising 
model depends on a parameter (3 which has to be estimated. More complex models depend on 
larger sets of parameters which we will denote by 9. 

The estimation of the energy parameters from an image is a difficult problem since it is 
not easy to apply standard estimation methods such as the maximum likelihood method. The 
likelihood function is given by 

p{I\6) = - log Z(0) -J2V C (I C ;0) (4.62) 

Unfortunately we do not have a closed form expression for the partition function nor an efficient 
way to compute it. The definition 

Z(0) = ^e-£<=€c^(ic;0) ( 4<63 ) 

cannot be applied in practice since it involves the sum of the non-normalized distribution for 
all possible image configurations. 
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4.9 Appendix - Fast Wavelet Transform 
Inner products of basis function 

Inner products of basis functions at consecutive levels < <^ jP , ^j-i, g >, < (j)j, 0 ^j-i,q > can be 
easily computed using the impulse responses ho, hi: 

< twtj-hq >= J ^ j x-p) ^~ l x - q)dx = V2 J <\>(2y + 2q - p) (j>{y)dy (4.64) 

< ^P^j-hQ >= J 2 '~ h <\>(Vx-p) ^~ x x - q)dx = V2 J 0(22/ + 2q - p) ^j(y)dy (4.65) 
where y — 2 J_1 x — q. Therefore, 

< <^\p> ^'-i,g >= h o(P ~ 2^) (4.66) 

< ^>,^j-i,g >= hip - 2q) (4.67) 

Fast wavelet transform: analysis 

The analysis step of the fast wavelet transform expresses the high resolution representation into 
a lower resolution representation. First, let us try to represent the basis functions of V}_i, Wj-i 
as a linear combination of the basis functions of Vj 

<t>j-i,k{x) = ^Cn^j^x) i/>j- lik (x) = ^dn^j^x) (4.68) 

n n 

where c n =< (j)j-i,k, 4>j,n >= h 0 (n - 2k), d n =< ipj-i,k, 4>j,n >= M n _ ( see (4.66,4.67)). 
The approximation and detail coefficients can be computed as follows 

Cj-i,k =< f,<l>j-i,k > d j-i,k =< f,*Pj-i,k > (4.69) 
Using (4.68) we obtain 

Cj-l, k =< /, 0j-l,k > = < /, Yj 2A 0<^> > (4.70) 

n 

n 

and 

=< f,ipj-i,k >=< /,$^ - 2fc)<^, n > (4.72) 

n 

= ^ /ii (n - 2&)c J> (4.73) 

n 

as wished. 
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Fast wavelet transform: synthesis 

Now we want to express (j)j n in terms of lower resolution basis functions 

= ^2°k<l>j-iA x ) + ^2 d k^j-i,k(x) (4.74) 



k k 



where c k =< ^jn^j-i^k >= h 0 (n - 2k), d k =< ^jn^j-i.k >= h 0 (n - 2k) Since, c J?n =< 
/, (j)j, n > we obtain 

c 3,n =< /, h °( n ~ 2k )<l ) 3-hk+ < /, 53 - 2fc)^_i, fc > (4.75) 

Cj, n = 53fto(n-2fc) </,^--i,fc >+53fti(n-2fc) <f^j- hk > (4.76) 
Cj, n = ^/i 0 (n - 2k)cj- 1:k + ^fti(n - 2k)dj- 1:k (4.77) 

as we wished. 
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Chapter 5 

Learning 



5.1 Introduction 

Learning methods play a central role in image and processing and vision and can be applied 
to a variety of problems ranging from texture analysis to image segmentation, object detection 
and tracking. 

Let us consider a couple of examples. Suppose we observe the motion of an object in an 
image. How can we predict the object position in the next frame and learn to do that in an 
automatic way? Suppose we observe an aerial image with different types of textures. How can 
we learn the main features of each region and find its boundaries? 

Learning can be done using deterministic or statistical methods. In both cases we want to 
predict a variable y given an observation x and try to do it by learning a function 

y = f(x,9) (5.1) 

which depends on an unknown parameter 9. 

We will address both approaches through a set of simple problems. 

5.2 Motion Prediction 

Suppose we observe the motion of a point object in an image (e.g., a circular motion) and want 
to predict the position of the object in the next frame. 

Let x(t) be the position of the point at frame t (t is an integer). We want to predict 
y(t) = x(t + 1). How can we learn this motion? 
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One strategy is the following. First we choose a model which we believe is rich enough to 
describe the motion 

y = f(x,9) (5.2) 

where 9 is a set of unknown variables. In the case of point motion we could assume for example 
that the predicted coordinates 2/1 2/2 (*) are related to the current coordinates xi(t),X2(t) by 



yxit) = aixi(t) + a 2 x 2 (t) + a 3 
y 2 (t) = a 4 x 1 (t) + a 5 x 2 (t) + a 6 



(5.3) 



with 9 — (ai . . . clq). Second, we need to learn the model parameters 9 from the available data. 
This can be done by collecting a set of input-output pairs {(x(t), y(t)), t — 1, . . . , N} and then 
adjust the unknown parameters 9 by fitting the model to the data using a sum of squared errors 
criterion 

£(0) = X>(t)-/(*(t),0)ll 2 (5.4) 



The parameters are estimated by minimizing this criterion 

9 — argmaxi£(0) 



(5.5) 



This is called off-line training since we first collect the data and fit the model later. We could 
also do both things simultaneously i.e., update the model every time we receive new information. 
In this case we would be doing do an on-line training. 

The third step consists of checking if the model is good enough. If it is not, we should try 
to improve it changing the model structure (5.2), the fitness criterion (5.4) or the optimization 
algorithm. 

This method is called the least squares method 1 and it is fully deterministic. 

Let us solve the prediction problem using the model 5.3. The model can be written as 



Vi (*) 




Xl(t) x 2 (t) 


1 0 


0 0 


2/2 (t) 




0 0 


0 xi{t) 


x 2 (t) 1 



(5.6) 



Hhe least squares method was proposed by Gauss in the 19th century, as a tool to predict the position of 
planets from experimental data 
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Figure 5.1: Motion prediction: original data (left), one step prediction (center) and prediction 
from the initial condition by recursive application of the one step predictor 



Writing the equations for the observed instants of time we obtain 



" J/l(l) " 




' *i(l) 


*2(1) 


1 


0 


0 


0 




1/2(1) 




0 


0 


0 


X1 (l) 


x 2 (l) 


1 




dl 


Vi(N) 




xi(N) 


x 2 (N) 


1 


0 


0 


0 






V2(N) 




0 


0 


0 


xi(N) 


x 2 (N) 


1 





(5.7) 



This can be written in a compact way as 

M9 = b 

The quadratic cost functional is given by 

E = (M<9 - b) T (M0 - b) 



(5.8) 



(5.9) 



and can be easily minimized by computing the derivative with respect to 6 and making it equal 
to zero. We obtain 

M T M9 = M T b (5.10) 

This is a system of 6 equations which can be easily solved. 

Figure 5.1 (left) shows 100 points of a circular trajectory with random disturbances. The 
motion model was then learned using this trajectory and the prediction estimates are shown at 
the center. To see how good is the motion, we generated the trajectory on the right recursively 
applying the predictor to predicted data. 



yi(t) = aiti(t - 1) + a 2 y 2 {t - 1) + a 3 
y 2 (t) = a 4 yi(t - 1) + a b y 2 (t - 1) + a 6 

This trajectory was generated without data. 



(5.11) 
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Least Squares Method 

Goal: predict a variable y from x 
Training data: {(x(t),y(t)),t = 1,...,7V} 

1. Choose the model: y — f(x,6), (9 unknown) 

2. Choose the cost function: E(9) = J2?=i \\y(t) - f(x(t), 0)\\ 2 

3. Estimate 9: 9 — argmax^ E(9) 

4. Validation: apply the model to new data and check if the desired accuracy was obtained. 

Table 5.1: Least squares method 
The least squares method is summarized in table 5.1. 

The solution is very simple because the cost function is quadratic and the model is linear 
with respect to the unknown parameters. We can use more flexible models with nonlinear terms 
in x (e.g., polynomials) 

p 

y = Y,a i f i (x ] 0) (5.12) 

1=1 

and still obtain a linear set of equations for 9. Things become more complicated when the 
model is nonlinear in 9 or the cost function is non quadratic. In these cases the optimization is 
often done using numeric recursive algorithms. 

5.3 Texture classification 

Let us consider a different problem. Suppose we have an image with two textures and want to 
automatically label each of them. We will make some simplifying hypothesis. We assume that 
the image / is a collection of N x N blocks each of them being filled by one of the textures. 

To solve this problem we first have to characterize each of the textures i.e., we have to learn 
their properties. For the sake of simplicity we assume that there are training images of the two 
textures and we can use use these images to learn the texture properties. Second we wish to 
estimate the label y G {0, 1} associated to each block of the test image. 

Figure 5.2 (left, center) shows the training images used to estimate the model. These textures 
were synthesized using the difference equation 

I(m, n) — ciio I (m — 1, ^) + aoi/(ra, n — 1) + aul(m — 1, n — 1) (5.13) 

with the following parameters: aoi = 0.4, aio = .l,&ii = .5 (texture 0), aoi = 0.1, aio = 
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Figure 5.2: Texture learning: two textures (left, center) and scatter diagram of texture features 
computed in 10 x 10 blocks 

.4, an = -5 (texture 1); w(m, n) is white (uncorrected) Gaussian noise with distribution N(0, 1). 
Figure 5.2 shows two images produced by this model. Both textures have zero mean and the 
same variance. These parameters are not able to discriminate the textures 

Both textures have directional structure. Therefore, we will adopt two directional features 
based on intensity differences along different directions 

Xl = E{(5| ?1 (m,n)} x 2 = E{S 2 12 (m,n)} (5.14) 

where Sij are intensity differences defined as follows: (%(ra, n) — /(m, n) — I(m — i,n — j). We 
compute a pair of features for each 10 x 10 block. The expected value is obtained averaging the 
squared differences in each block. Figure 5.2 (right) shows 100 feature points for each texture 
computed in 10 x 10 blocks. This scatter plot shows that we can discriminate both textures 
with few errors since the overlap of the two clouds is small. 

The next question is: how do we classify texture blocks using this information? If we knew 
the probability distribution of the features associated to each class the answer would be simple. 
The optimal decision is based on the Bayes law 

P(y\x) = ap(x\y)P(y) (5.15) 

where P(y) is the a prior probability of the label y, p(y\x) is the a posteriori probability of 
the label y, p{x\y) is the data distribution and a is a constant. Note that y takes two possible 
values: 0, 1. The classifier should choose the the most probable label i.e., the label y which has 
the highest probability P(y\x). This classifier is called the Bayes classifier and it is optimal in 
the sense that it has the smallest probability of error. 
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Figure 5.3: Texture classification: input image (left) true class (center) estimated class (right) 

In many problems the probability distributions are unknown and all we know is a set of 
input-output pairs (xi^yi). In this case we have to learn the classifier and the story is much more 
complex. In fact there is a whole area in statistical learning which deals with this problem known 
as supervised learning. Some popular methods to classify data from a training set are: Fisher 
discriminant, the nearest neighbor classifier, classification trees, and support vector machines 
[3]. 

We will adopt a different approach. We will estimate the data distributions p(x\y) and use 
the Bayes classifier. The data distribution are Gaussian N(fi yj R y ) J y G {0,1} where fi y ,R y 
are the mean vector and covariance matrix which can be estimated from the training set by 
standard formulas 

1 N 1 N 

V=j^Ys Xi R = j^Y^( Xi ~ V)( x i~ V) T (5.16) 

i=l i=l 

Using the data points in Figure 5.2 (right), we obtained 



3.57 




1.80 




1.103 0.068 




0.126 0.147 








Ri = 








1.73 




3.40 




0.068 0.127 




0.147 1.852 



(5.17) 



The most probable labels can now be computed using the Bayes law (5.15). The normal 
probability density function is 



7V(x;/i,R) 



-i(x-^) T R- 1 (x- A x) 



(5.18) 



(27r) d / 2 |R| 

where d = 2 is the dimension of x and |R| is the determinant of the covariance matrix. Figure 
5.4 shows the test image made of two types of textures (left), the true label of the image blocks 
(center) and the estimated label (right). The estimated labels are identical to the true ones. No 
labeling error was made since the textures can be easily discriminated. 

The Bayes classifier (known distributions) is summarized in table 5.3. The strategy followed 
in this problem (we used the expression of the Bayes classifier replacing the unknown parameters 
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Bayes Classifier 

Goal: assign a label y to an observation x 
Hypothesis: p(x\y), P(x) are known 

1. Compute the a posteriori probabilities: P(y\x) — ap{x\y)P{x) 

2. Choose the most probable label y 

Table 5.2: Bayes classifier 

by estimates) is suboptimal but it is a very common practice. This method is known as the 
plug-in Bayes classifier 2 . 

In the classification problem solved in this section we made some simplifying hypotheses. 
We will try to relax some of them in the next section. 

5.4 Parameter estimation 

We had to estimate the parameters of a multivariate Gaussian distribution from data in the 
previous examples. Fortunately, there are well known formulas for the mean and covariance 
matrix in this case and we did not bother to derive them. How should we proceed when we do 
not know expressions for the parameters estimates? 

There are several well known methods to solve this problem. We will briefly review two of 
the most popular: the maximum likelihood methods proposed by Fisher and the MAP method. 

Maximum Likelihood Method 

Suppose we want to estimate a probability density function p{x\9) of a random variable x where 
9 is an unknown parameter. In general x and 9 are vectors. Let X — {xi,i — 1, . . . ,n} be 
a set independent realizations of x known as training set. We want to estimate the unknown 
parameter 9 using the X. 

The maximum likelihood (ML) method chooses the 9 which makes the observations more 
probable i.e., which maximizes p(X\9). The estimate is therefore defined by 

9 = arg max L(9) (5.19) 

2 the plug-in Bayes classifier is suboptimal because it neglects the uncertainty of the parameter estimates 
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where 

n 

L(e)=p(X\8) = l[p(x i \8) (5.20) 

is known as the Likelihood function. 

Sometimes we use the log likelihood function instead 

0 = axgmaxZ(0) (5.21) 

n 

1(0) = \ogp(X\0) = J2^gp(xi\0) (5.22) 

which is often simpler. 



Example 

Suppose we wish to estimate the parameter s of a Rayleigh distribution 



3 



P {x) = ^e~%* (5.23) 



from a set of independent realization X — {xi}. 
The log likelihood function is given by 



/ ( s ) = Ei e "^ ( 5 - 24 ) 
i=i s 

Therefore, 

l(s) = C-2n\ogs- ^S 2 (5.25) 

where 

S 2 = l±xl (5.26) 

i=l 

Computing the derivative of l(s) and making it equal to zero we obtain 

-j +n-^S 2 = 0 (5.27) 

s 2 = ^S 2 (5.28) 

This is the ML estimate of s. 



3 the Rayleigh distribution is used to describe the intensity of ultrasound images. In this case the parameter 
s varies along the image. 
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Maximum a Posteriori Method 

The maximum a posteriori (MAP) method is a Bayesian method. It assumes that 9 is itself a 
random variable and we know its distribution p{9) also called the a priori distribution or simply 
prior. 

After observing the data X we can update the a priori distribution of 9 using the Bayes law 
and obtain 

p(0\X) = ap(X\9)p(9) (5.29) 

where a — p(X). This is the so-called a posteriori distribution. Note that there is a full 
equivalence with what we did before for the Bayes classifier. Everything is the same except 
that we now wish to estimate a parameter (integer or a vector) instead of a class label. 

The MAP method chooses the most probable value of 9 using the a posteriori distribution 

9 = argmaxp(<9|X) (5.30) 

9 = arg max p(X 1 9)p(9) (5.31) 

We note the similarity with the ML method. The MAP method maximizes the ML function 
multiplied by the prior. The prior is the new term which accounts for the existing information 
about the parameter 9 before the observations. 

This is also a sound framework to perform data fusion if we have multiple sources of infor- 
mation or information obtained at different instants of time. The estimate available at time 
t — 1 can be updated with new information X 1 by 

9 l = argmaxp(X*|<9)^-i(<9) (5.32) 

9 

where 

Pt - 1 (0)=p(0\X\...,X t - 1 ) (5.33) 
is the distribution of 9 with all the information available until time t — 1 and updated by 

Pt(9) = atPiX^pt-iV) (5.34) 



Example Suppose we want to estimate the mean of a normal distribution N(fi J a 2 ) with 
known a known that the prior distribution of fi is normal N(/jlq, ctq). 
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The log likelihood is 

n -j 

l(f,) = logH^^e-^-rf (5.35) 



1=1 v 



1 " 

/(/*) = -nlogA/W-- T ^(a ;i -/x) 2 (5.36) 



2 ^ 2 i=1 

The logarithm of the prior is, 



logp(p) = — log27r<7 2 - — {p - mo) 2 (5.37) 



ilog2™ 0 2 -^ 
The function to be maximized is 



1 n 1 
+ logp(/i) = C - ^2 " ^ " ^2 " ^? (5.38) 



Computing the derivative with respect to fi we obtain, 

1 " 1 

— J2(fi-xi) + — ( M -Mo) = 0 (5.39) 
a <=i a o 

— + — U= — + % 5.40 
a 2 a 2 / a 2 a 2 

S + ^ 

where x is the average value of x in the training set. The estimate obtained with the MAP 
method depends on the data and on the prior. The influence of the prior is negligeable if n is 
large. However if we have few data points the prior plays an important role. 



5.5 Unsupervised Texture Classification 

Suppose we had a textured image to learn the classifier but we did not know the true label of 
each block. The Data used to learn the classifier only has the input x but does not include the 
desired output (label) y. The training set is therefore given by {xi}. How do we estimate the 
relationship between x and y if we do not have any examples of y. 

This seems an impossible problem since we know nothing about the input-output (x,y) 
relationship. 
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It is amazing but we can still solve this problem in an efficient way and achieve good results 
if the two clouds do not overlap too mutch. Suppose we have the same training data as before 
but without the labels i.e., we have the data points shown in Fig. 5.2 (right) but they do not 
have the color labels since the label (color) is unknown. We can still separate the two clouds 
and estimate the mean vector and covariance. This seems intuitive if we look at the figure and 
look for ellipsoidal clouds of points. However, how do we do this in an automatic way? 

The probability distribution is a mixture of two Gaussians 

p(y) = c 0 N(x; /i 0 , Ro) + ci7V(x; ^ , R ± ) (5.42) 
N(x;^R) = _jl_ e -i(«-/') T «- 1 («-/') (5.43) 

where Co, ci(co + c\ — 1) are the mixing probabilities and fii^Ri^i — 0,1 are the parameters of 
the Gaussians. How can we estimate the mixture parameters? 

This problem is solved by the Expectation-Maximization (EM) method. This method allows 
us to fit the data set {xi} by two Gaussian distributions. The idea is simple. First we initialize 
the mean and covariance of the two Gaussians (/ii , \i<i , R\ , R2) . Then we compute the probability 
of both labels for each data point i.e., we compute the weights 

Wij = P(yi\xi) = ap(xi\y = j)P(y = j) (5.44) 

The weights sum one (wio + wn = 1) and they have an intuitive meaning: is the degree of 
membership of Xi to label j. After this step, we update the Gaussian parameters using all the 
data points wighted by their associated probabilities as follows 

= Eili WijXi R = Eili w ij( x i - v)( x i - n) t / 5<45) 
E i= i w ij E i= i w ij 

Let us apply the EM method to the data points of Figure 5.2 (right) without labels (colors). 
The mean vectors \i\ , /12 were initialized by randomly selecting two data vectors and the co- 
variance matrices Ri,R2 were made equal to the identity matrix. Figure 5.4 shows the results 
obtained at iterations 0 (initialization), 1,5 and 20. It shows level curves of the probability 
distribution which are ellipses and the weights associated to the label 0. Initially the weights 
are close to 0.5 for most of the data points because the covariance matrices are very large and 
all data vectors can be represented by both Gaussians. As the estimates improve the weight 
vectors become closer either to 0 or 1 as expected. 

After estimating the probabilistic model using the EM method we applied this model to 
classify the texture blocks of the test image 5.4. We have obtained a perfect classification as 
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Figure 5.4: Mixture estimation with the EM method: data points, level curves and degree of 
membership woj at iterations 0 (random initialization), 1, 5, 20. 

before. The EM method for a mixture of Gaussians is summarized in table 5.3. The method 
can be used in under more general conditions [3]. 



5.6 Dealing with Arbitrary Shapes 

We discussed how two label different textures by learning a function y — f(x;0) when we have 
a labeled training set (supervised learning) and when we do not known the correct labels in 
the training set (unsupervised learning). We have assumed that the image consisted of non 
overlapping blocks of homogeneous texture. 

What should we do when this hypothesis is not valid and the labeled regions have arbitrary 
shapes? This is a difficult step. First feature extraction becomes more difficult since it depends 
on the labeling step. We can still compute the mean of the intensity differences in each region 
but the regions shapes depend of the labeling. 

Second we need a label for each pixel instead of assigning a common label to all the pixels 
of a block. Region shape influences feature extraction and vice versa. 

We will discuss this type of problem latter using random space models (Markov random 
fields). 
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EM method 

Goal: estimate the parameters of a Gaussian mixture p(x) — Y^JjLi c jN{Hj,Rj) 
Data: {xi} 

initial 
iterate: 

E step: compute the label probabilities for each data point 

Wij — aiN(x; [ij,Rj)cj j — 1, . . . , m (5.46) 

M step: update the mixture parameters (cj^fij^Rj) 

c, = ^f>, H = ^f^ Rj = YL^-^-Hf (5 . 47) 

iV <=i E<=i w a Ei=i w a 

end 



Table 5.3: EM Algorithm: estimation of mixtures of Gaussians 
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Part II 

Image Processing 
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Chapter 6 

Image Enhancement and 
Restoration 

6.1 Introduction 

Images are often degraded during the acquisition process (camera out of focus, moving camera, 
lost image pixels, noise). The goal of image enhancement and restoration is to recover from 
these errors as much as possible. 

To perform this operation we need to make hypothesis about the image properties and 
about the type of degradation that occurred. These assumptions may be specified by the user 
or learned from the data. For example, if the observed image is corrupted by noise we can either 
specify the noise statistics (e.g., white Gaussian noise) or learn the noise model from original 
and distorted images. 

This topic is an excellent introduction to model based image processing since it combines 
many important issues: image models, acquisition models and inference techniques. The meth- 
ods used in this section are useful to address other image processing problems such as image 
alignment and super resolution. 
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Figure 6.1: Degradation with white Gaussian noise 

6.2 Inverse problems 

Image restoration is an example of a wider class of problems called inverse problems. These 
problems aim to invert the transformation 

y = Ax (6.1) 

where y is the observed variable, x is the unknown variable and A is a non invertible operator 
or an ill-conditioned operator. The above equation does not have a unique solution. 
Let us consider the equation 

y = x 1 -\-x 2 (6.2) 

where xi,x 2 are two unknowns. It is not easy to determine them from their sum y. 
This chapter gives several methods to address this type of problems. 



6.3 Degradation Model 

We will assume that there is an original image (unknown) which is distorted in some way. For 
example, let us assume that the original image /(x) is corrupted by white Gaussian noise. The 
observed image is therefore 

J(x) = J(x) + W(x) (6.3) 

where W(x) is a noise image i.e., £{W(x)W(x')} = cr 2 £(x, x'), W(x) ~ J\f(0,a 2 ). Figure 6.1 
shows an example of this degradation process. 

The acquisition process is often more complex involving two types of degradation: a deter- 
ministic degradation (e.g., burr) and a random degradation. In this case we can often assume 
that 

J = T(J) + W (6.4) 
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where T is a deterministic operator and W is the observation noise. In some problems we may 
assume that T is a linear operator. In this case, J is the convolution of the original image / 
with an impulse response: 

T(/)(x) = ^ /i(x, x / )/(x / ) space- varying (6.5) 

x' 

or 

T(/)(x) = ^ /i(x — x / )/(x / ) space-invariant (6.6) 

x' 

For example, out of focus picture and motion blur can be modeled in this way using appropriate 
degradation filters. In the first case the impulse response is isotropic while in the second the 
impulse response depends on the motion direction. 

It is often useful to represent the original and observed images with vector notation. Let I, J 
be the vectors obtained by concatenation of the columns if the original and degraded images; 
the dimension of these vectors is n x 1 where n is the number of pixels of the image 1 . 

The linear degradation process can now be stated as 

J = MI + W (6.7) 

where M is a degradation matrix and W a random matrix. This is a very useful notation. 
However, it should not be used to simulate image degradation in a computer since the dimension 
on M is very high (e.g., 10 6 x 10 6 ). 



6.4 The naive approach 

To estimate I from J we must be able to separate MI from W in (6.7). This is only possible 
if we make some hypotheses about both terms. 

Let us ignore this statement for the moment and try to apply a classic estimation method 
(maximum likelihood) in a blind way. 

Assuming that the noise image W is white and Gaussian the log likelihood function defined 

by 

1(1) = logp(J|/) (6.8) 

can be easily computed 

1(1) = lo g ;Q^(x)|/(x)) = 5>gp(J(x)|J(x)) (6.9) 

X X 

Hhe symbol I is being used with two meanings: original image and identity matrix; the meaning should be 
clear from the context 
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Since the noise is Gaussian p(J(x)\I(x) — N(J(x); I(x),a 2 )) 

1(1) = -| log(27ra 2 ) - Y, ^2 ( J W - / W) 2 ( 6 - 10 ) 

X 

Function / is maximized if /(x) = This means that the ML estimate of / is the 

observed image J. Discouraging. According to this method we should not modify the observed 
image. What is missing? The maximum likelihood model makes no assumption about one of 
the signals: the original image. We need to use a model for the image as well. 

In image restoration problems there is often another difficulty . There is also loss of infor- 
mation when we apply the operator T. The operator is non invertible is many problems. 



Example 

Let us consider the following problem: I = (h,h) is an unknown vector and we only observe 
the sum of both components (J = 1\ + I2). There is also loss of information in this case. How 
can we estimate I from J? 

We must make some assumptions about the solution. 



6.5 Regularization Methods 

The previous difficulties can be overcome by assuming that the reconstructed signal / has to 
verify additional restrictions. This raises an important question: how can we characterize the 
class of admissible images. The answer to this question has a direct influence on the solution. 

In this section we will assume that the image / is smooth. We wish to find an image / such 
that T(I) ~ J i.e., the (squared) norm of the reconstruction error 

\\J-T(I)\\ 2 (6.11) 

is small. However this is not enough to obtain an unique solution and we wish to impose an 
additional constraint. We wish to find the smoother image among all the images with small 
reconstruction error. The smoothness criterion is defined as 

\\D(I)\\ 2 (6.12) 
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where D is a differential operator, e.g., an operator which computes difference between pairs of 
neighboring pixels 2 . 

These two different goals can be included in a single cost functional 

C=\\J-T(I)\\ 2 + a\\D(I)\\ 2 (6.13) 

with two terms: a data term which measures the fitness of the model to the observed data and 
a regularization term; a controls the relative importance of both criteria. The minimization of 
C leads to a reconstructed image which is compatible with the observations and it is smooth. 
The methods based on the minimization of (6.13) are called regularization methods. 

These principles can be applied to the denoising problem using the Euclidean norm. In this 
case, 

C = £(J(x)-/(x)) 2 +a Y, '(y)) 2 (6-14) 

x (x,y)€V 

where V the set of all neighboring pixels. 

To minimize C, all the partial derivatives with respect to J(x) must be zero 

5C 



<5/(x) 

Therefore, 



= 0 (6.15) 



-2( J(x) - /(x)) + a 2(/(x) - I(y)) = 0 (6.16) 
yev x 

where V x is the set of neighbors of pixel x. Therefore, we obtained a system of linear equations 

(1 + aiV x )/(x) - a Hy) = ^(x) (6.17) 

where 7V X is the number of neighbors of pixel x. The previous equations can be recursively 
solved. 

where /(x) is the average intensity of all the neighbors of pixel x (excluding x), computed using 
the last estimate of the image. Pixel update may be performed using two different strategies. 
We may update each pixel as soon as we have a new estimate for it. In this case we are solving 
the system of equations using the Gauss-Seidel method. Instead, we may compute all pixel 
estimates first and update all the pixels simultaneously. In this case we are using the Jacobi 
method. 



2 we note that the two norms are defined in different spaces and can be independently chosen. 
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Matrix notation is very useful to describe the regularization method with quadratic cost 
even if some of the matrices can not be stored on a computer due to their huge dimensions. Let 
I, J be two column vectors obtained by stacking the columns of the corresponding images. The 
degradation model can be written as J = I + W where W is a noise vector. The cost functional 
can now be written as 

C = ||J-I|| 2 + a||DI|| 2 (6.19) 

C = (J - I) T (J - I) + cd T D T DI (6.20) 

where ||.|| stands for the Euclidean norm and D is a sparse matrix which computes the differences 
between pairs of neighboring pixels. Computing the derivative of C with respect to vector I we 
obtain 3 

-(J - I) + c*D T DI = 0 . (6.21) 

and 

(I + cdD T D)I = J (6.22) 

This is a linear system of equations written with matrix notation. 

Although this equations is simple it is not easy to invert matrix I + aD T D since it has a 
huge dimension nxn where n is the number of pixels in the image. In fact we cannot even store 
it in a computer. In practice, recursive algorithms are used to solve the system of equations 
such as the ones described before. 

Figure 6.2 (top-right) shows the cameraman figure corrupted by white noise with SNR=13dB. 
It also shows the reconstruction results obtained with the regularization method using a — 1 
(SNR=17.0dB) and a — 8(SNR=18.6dB). The first choice of a has corresponds to a mild reg- 
ularization effect. The input noise is still clearly visible at the output. The second choice for a 
corresponds to a strong regularization which attenuates the noise but it also blurs the image. 
The final result looks as if the image was out of focus. This is the trade off involved in the 
choice of a. 

6.6 Edge preserving methods 

The previous method attenuates the noise by performing a low pass filtering of the signal. 
Therefore, fast transitions are filtered and replaced by smooth transitions. This creates a burr 
effect which is perceptually annoying since it destroys the edge information. 

3 see the derivative rules in Appendix 
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Figure 6.3: Cost functions: x 2 , \x\ and x 2 /(l + x 2 ) (x and y scales are different) 

This raised an important question in image reconstruction: can we attenuate the noise 
without destroying the edges? 

The bad performance of the previous method at the transitions is a consequence of the way 
transitions are penalized by the cost function C. The squared Euclidean distance ||/(x) — J(y)|| 2 
rapidly grows when there is a sudden change of intensity in image /. It would be better if we 
the transition cost grew slower (see figure 6.3) in such a way that small changes due to noise 
would still be penalized but very large changes due to transitions would not receive a very high 
penalty. 

One way to reduce the cost of abrupt changes consists of using the l± norm in the regular- 
ization term 

c = 5>(x)-/(x)) 2 + a l'( x W(y)l ( 6 - 23 ) 

x (x,y)GV 

This norm grows slower than the I2 norm and it is less sensitive to abrupt changes. However, 
the minimization of C becomes a nonlinear problem. 

This difficulty can be partially overcome by converting this problem into a sequence of 
quadratic problems which can be solved by linear methods. Equation (6.23) can be written as 



where 



C = £(J(x)-J(x)) 2 +a Yl Mx,y)(/(x)-/(y)) 2 

* (x,y)ev 



w(x,y) = 



(6.24) 



(6.25) 



|/(x)-/(y)|+e 

is a weight and e is a small constant used to prevent the weights from tending to infinity. The 
minimization of (6.24) is similar to the minimization of (6.23) if we consider the dependence of 
the weights w(x, y) on the image I(x) — /(y). However, if the weights are kept constant in each 
iteration (simplifying hypothesis) the problem becomes mutch simpler since the cost functional 
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Figure 6.4: Denoising with the regularization method: (a) original and observed signals; restored 
signals with (b) squared Euclidean norm (SNR=17.4 dB) and with (c) the l\ norm (SNR=21.7 
dB) 



becomes quadratic and can be solved as follows 



[i + a ™( x > y)] = + a ^ ™( x > y)Hy) 



/(x) 



j(x) + a E yG y x ^( x > y) J (y) 



(6.26) 



(6.27) 



1 + aE yG v x ^( x >y) 

Since the weights depend on the estimate / they must be updated at the end of each iteration 
and the cost functional must be minimized again. The process usually converges after a few 
iterations. 

Figure 6.4 illustrates the performance of the regularization method at transitions using the 
norms I2 and l\ . The best results at transitions are achieved with the norm li while the I2 norm 
achieves slightly better results at smooth regions. It is also observed that the l± norm tends to 
produce a piecewise constant estimate i.e., it is very well adapted to that class of signals. The 
output of this algorithm in the cameraman example is shown in Figure 6.5. 



6.7 Bayesian Methods 

The Bayesian methods formulate the reconstruction problem in a different way but they lead 
to similar algorithms as we shall see in this section. 

Bayesian methods assume that the unknown image / is a realization of a random variable 
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Figure 6.5: Denoising with the norm L\ (left, SNR=19.3dB) and using wavelets with soft 
thresholding (right, SNR=18.5dB) 

with a known distribution p(I) called the prior. The prior contains all the available information 
about the statistical properties of the image e.g., that it is a smooth signal in most of the points 
and may exhibit sudden changes at a small subset of points. 

Bayesian methods compute the a posteriori distribution of I after observing the noisy image 
J. The a posteriori distribution is obtained using the Bayes law 



p(I\J) = kp(J\I)p(I) 



(6.28) 



where p(J|I) is called the sensor model and describes the degradation process, p(I) is the prior 
and k is a normalization constant. 

To estimate I we can compute the maximum of the a posteriori distribution p(I|J). This is 
the most probable value for the unknown image and the method is known as the maximum a 
posteriori (MAP) method 

f = argmaxp(I|J) (6.29) 

The MAP method requires the specification of p(J|I) (sensor model) and p(I) prior distri- 
bution for the image. The method may lead to simple or complex a posteriori distributions 
depending on the sensor model and prior. A simple case is obtained if the image is corrupted by 
white Gaussian noise W ~ Af(0, cr 2 I), and the prior is Gaussian as well I ~ Af(0, (aD T D) _1 ). 
In this case 

Pim = 6XP( - 2^ (J - I)T(J - I)} (6 - 30) 
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eM~^ T D T DI) 



(6.31) 



The image estimate I is obtained maximizing 



p(J|I)p(I) = Kexp(-^(J - I) T (J - I) - |l T D T DI) 



(6.32) 



This is equivalent to the minimization of 



C=(J-I) T (J-I)-a 2 aI T D T DI 



(6.33) 



which is the same as the regularization cost functional (6.20) with Euclidean norms. Both 
methods become identical if the regularization term is based on the Euclidean norm and the 
prior is Gaussian prior iV(0, (aDD T ) _1 ). The prior plays the role of the regularization term. 

6.7.1 Denoising with wavelets 

Explicar o soft e hard thresholding. 



6.8 Moisaicing and Superresolution 
Exercises 

• Extend the regularization method for the linear model J = HI + W, where H is a square 
image. 
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Chapter 7 

Image Segmentation 



7.1 Introduction 

Image segmentation is a basic image operation. It amounts to finding regions associated with 
the objects of interest, in the image domain. Our visual system performs this task in a natural 
way. We have no difficulty to distinguish all sorts of objects when we look at a scene. However, 
segmentation is one of the most challenging operations in image analysis. 

Figure 7.1 shows some examples from different application areas. Segmentation is used in 
document image analysis to separate written text from the background, allowing the application 
of recognition modules. In remote sensing, image segmentation is used for a variety of problems 
e.g., to classify the terrain cover as forest, agriculture, urban, water. Image segmentation is 
also a key operation in medical applications. It is used to extract the boundary of organs and 
anatomic structures in images obtained by different types of modalities (CAT, MRI, US). A 
fourth example is drawn from the video surveillance area. Modern video surveillance systems 
must be able to segment moving objects in an image (e.g., vehicles in highways) in order to 
classify them as well as their activities (e.g., classify dangerous maneuvers is a highway). 

Image segmentation is a difficult task. It is not easy to include the a priori information we 
have about the object (shape, color) in an objective criterion. Most algorithms are still far from 
that perspective and are based on two simple ideas: 

• region homogeneity: the image should be split into homogeneous regions with similar 
intensity, color, texture or motion 

• transitions: transitions (e.g., edge points) detected in the images are used to estimate 
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Figure 7.1: Segmentation problems 



the object boundaries. 

The first type of methods estimate regions using an homogeneity criterion while the second ones 
estimate the objects boundaries. There is an underlying assumption that objects are associated 
homogeneous regions with abrupt transitions at their boundary. This is not true in many cases. 
Therefore the performance of most methods is still limited. 

Prior information about the scene plays a very important role. If we are segmenting cows 
in a green field we should take this information into account. Therefore segmentation is not 
a pure bottom up problem. It should take into account the available information about the 
objects present in the scene (top-down). This is a challenging problem because it is not easy 
to represent the objects shape, color and texture information in useful way. How can we teach 
a computer to segment people in a mob? 

The next sections formulate the segmentation problem and present some of the available 
techniques from simple thresholding methods to Markov random fields or level set approaches. 

7.2 Problem Formulation 

Given an image u we wish to split the image domain into n regions Sk,k = 0, . . . ,L — 1 where So 
is the background image and the other regions are associated with different objects of interest 
(Figure 7.2). We will assume that {Sk} is a partition of the image domain i.e., all Sk are disjoint 
and their union is the whole domain. 
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Figure 7.2: Homogeneous regions 



This operation can be performed by assigning a label 



l k £ {0,l,...,n-l} 



(7.1) 



to each pixel (k,u(k)). The label sequence / = (Zq, • • • , Il-i) is equivalent to the image partition 
into disjoint regions {X^}. 

The human visual system performs visual segmentation using many types of cues and a 
great amount of knowledge about the visual appearance of different types of objects in different 
poses. This is a formidable task which can not be performed by current image analysis systems. 

Most of the approaches described bellow are based on some kind of model or homogeinity 
criterion which is assumed to be valid to represent the image in each region Sk . For each specific 
application we must be able to choose the most appropriate method. 



Consider the image u shown in Figure 7.1a. The text is darker than the background. One way 
to separate the text from the background consists of comparing the intensity of each pixel with 
a threshold T and assigning the label 1 if the intensity is smaller than the threshold. Therefore, 



Figure 7.3 shows the binary segmentation obtained with this method. This technique is called 
thresholding and it is used when we need a fast binary segmentation under simple conditions 
(the two classes are separable in the feature domain). 



7.3 Thresholding 




J(k) < T 



otherwise 



(7.2) 



105 




Figure 7.3: Binary segmentation obtained by thresholding 



7.4 Pixel Classification 

We can extend these ideas to L-ary decision problems. Suppose we wish to separate an image 
/ into L regions and we can assume that the pixels inside each region Si are independent 
realizations of a random variable with probability distribution p(I(k)\0i) where Oi characterizes 
the color or intensity of the i-th region (see Figure 7.2). 

Since the pixels are independent we can classify each of them in an independent way by 
choosing the most probable label Ik which maximizes the class probability 



where Pi is the a priori probability of the i — th class and a — l/p(/(k)) is a normalization 
constant which does not influence the decision. 

Instead, we could assign a cost to each label as follows 



and choose the label with smallest cost. 

The key point here is how to obtain the models p(I\0i) associated to the different objects 
(regions). A similar problem was discussed in chapter 5. 

We can distinguish two approaches. In the first approach the user is asked to identify the 
regions Si associated to each class of objects. In this case we have a supervised learning problem. 
The classified images provided by the user allow us to estimate the unknown parameters Oi using 
standard estimation algorithms e.g., the ML method 



P(i\IQs))=ap(IQs)\8 i )P i 



(7.3) 



£*(0 = -logP(/(k)|e«)-log^ 



(7.4) 



Oi = argmaxp(/ i |6» i ) 



(7.5) 



0 i — arg max 



II Pimm 



(7.6) 



k€Xi 
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where Ii — {/(k),k G Xj}. Suppose we want to classify the land cover using satellite images. 
First we would select patches from all the classes (e.g., water, agriculture, urban areas) to be 
separated. Then we would use learning methods to estimate the statistical models and classify 
the other pixels. This procedure can be applied to many problems (e.g., motion detection, 
pedestrian detection, face detection). 

When we do not have classified data the learning problem becomes much more difficult 
since the image contains L objects but we do not know where they are. This is an unsupervised 
learning problem which can be solved using clustering algorithms [3]. The distribution of the 
data is a mixture 

L-1 

p(I(k)\c,0) = J2c iP (I(k)\0 i ) (7.7) 

i=0 

where c = (co, . . . , cl — 1) are the mixture coefficients and 9 — (Oo, . . . , #l-i) are the parameters 
of the mixture modes. The mixture coefficients are the a priori probabilities of the classes. 
Therefore c& > 0 and 

L-1 

E c * = 1 ( 7 - 8 ) 

k=0 

The mixture mode p(/(k)|0i) is the distributions of the data associated to the i-th class. The 
estimation of c, 0 is much more difficult than in the supervised case because we do not know 
what is the class of each observation /(k). 

The Expectation-Maximization (EM) algorithm is a very useful technique to estimate the 
mixture parameters. The EM method is an iterative algorithm with 2 steps in each iteration. 
In the E-step, all the data points are classified in one of the available classes. This is done using 
a soft classifier i.e., by computing the probability of each class 

wj U = P(lk = i\I(k),0)) (7.9) 

Wk i = ap(I(k)\0 i )c i (7.10) 

using the current parameters c. The parameters of the mixture are updated using available 
weights computer by the EM method (see Chapter 5). 

Background subtraction 

Figure ? shows the image obtained with a surveillance camera and the evolution of a the 
intensity of a given pixel during a time interval. Most of the time the intensity is close to an 
average intensity /(k) plus some additional noise. Therefore, 

J(M) = /(k)+w(M) (7.11) 



107 



Figure 7.4: Motion detection by background subtraction [2]. Active regions are represented by 
the corresponding bounding boxes 

where w(k, t) is a white Gaussian noise with zero mean and variance cr(k) 2 . Therefore, w(k, t) ~ 
iV(0, cr(k) 2 ). The mean and variance of each pixel can be computed from a video sequence using 
robust estimation methods. The image / is often called the background image. 

Using this model, we can detection motion by comparing the absolute error |/(k,t) — /(k)| 
with a known threshold. Motion is detected in pixel k if 

|I(M)-/(k)|>»7 (7.12) 

This method is known as background subtraction method and it is often used in real time 
surveillance operations. 

This method is simple and fast. However it has two major drawbacks. It can not cope with 
intermittent noise (e.g., moving trees) and it can generate ghost active regions is a background 
object (e.g., a stopped car) starts to move. 

More sophisticated models are often used to cope with the first problem. For example, the 
pixel intensities are often modeled by mixtures of Gaussians which are able to represent multi 
modal distributions. This is very important to represent more complex patterns of motion as 
the ones generated by moving trees. The mixture parameters are often updated during the 
operation of the system to deal with changes of illumination. 

The independent classification of the pixels often leads to small groups of active pixels 
which can not be interpreted as objects. The binary image produced by the motion detector 
is often processed by morphological filters to remove small regions and fill the holes. This 
operation is performed in a suboptimal way since the post-processing criteria are different from 
the segmentation criterion. This difficulty can be solved by performing a joint classification of 
the image pixels. 
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7.5 Joint Pixel Classification 



The section addresses the joint classification of all the pixels: can we enforce labeling coherence 
to avoid spurious labels? The answer is yes and it leads to powerful (but time consuming) 
segmentation techniques. The basic idea is to include some prior information about label 
distribution. 

A simple way to encode this information is by saying that we want a small number of label 
changes among neighboring pixels. This can be done by defining a prior distribution P(x;a) 
for the label sequence / = (7o, • • • , l n -i) [pc represents a parameter to be estimated). This can 
be done using Markov random fields (see Chapter 4) e.g., the Ising model. 

Assuming a prior distribution, image segmentation can be performed by computing the most 
probable labeling i.e., by maximizing 

P(l\u) = cp(u\l) P(l\a) (7.13) 

This approach allows important improvements is most problems. However, the optimization 
of P(l\u) is usually a difficult combinatorial problem. 

The segmentation of the image can also be formulated in a deterministic framework. Let us 
suppose that there is a classification cost E k {l k ) when we assign a label l k to the k — th pixel 
and let Ei nt (l) be the cost of a labeling configuration x. For example, labeling configurations 
with a larger number of transitions will have larger costs. 

The energy to be minimized has two terms, a data term and a regularization term 

n-l 

E(l) = Y,E k (l k ) + E int (x) (7.14) 

k=0 

This approach is equivalent to the the most probable labeling (7.13) if we define 

E(l) = -logP(l\u) (7.15) 

In this case 

E(l) = - log P(u\l) - log P(l\a) (7.16) 

n-l 

E(l) = - log P(u k \l k ) - log P(l\a) (7.17) 

k=0 

Comparing with (7.13) we get an expression for the data energy and for the internal energy 

E k {l k ) = -\ogP{u k \l k ) (7.18) 
E int (x) = -log P(l\a) (7.19) 
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Figure 7.5: Watershed Transform: a) two types of drop locations b) catchment basins and 
watershed lines 

The two approaches (deterministic and probabilistic) are therefore equivalent in this case. 

One possible choice for the internal energy is to make it proportional to the number of label 
changes among all neighboring pixels 



where J\f is the set of all pairs of neighbors and 5 is the Kronecker symbol (£(i, j) = 1 if i = j 
and 8(i,j) — 0 otherwise). 

The optimization of E is a non-trivial task since it is an combinatorial problem with a very 
large number of variables. The ICM algorithm and the Gibbs sampler are popular choices for 
many years. However they have been outperformed by min cut / max flow algorithms which are 
able to obtain the optimal segmentation in polynomial (quasi-linear) time. In the binary case, 
these algorithms achieve optimal solutions for the segmentation task with a small computation 
time and they have already been used in real time surveillance applications. 

7.6 Watershed Transform 

The watershed transform (WT) is inspired in the ways water appears and spreads in nature. 
Suppose we interpret the image as a topographic map. The image intensity I{x) stands for the 
elevation at the point of coordinates x and suppose we place a drop of water on each point of 
the surface. Most of the drops will slide along the surface until they reach a local minimum (see 
Fig. 7.5). Few drops will not move since they are located on a local maximum. These drops 
will converge to more than one minima if they are slightly disturbed. 

We denote by watershed or catchment basin the attraction region of each minimizer x* of the 




(7.20) 



p,qEAf 
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Figure 7.6: Watershed segmentation: a) original image, b) gradient intensity, c) watershed 
segmentation 

image /. If we place a drop of water at each point of a catchment basis they will converge to the 
same minimizer. The points which become attracted to multiple minima if a small disturbance 
is applied to the drop location belong to watershed lines. 

The set of all catchment basins defines a segmentation of the image. Each region (object) 
is associated to a local minimum of the image intensity. 

The WT is a useful tool if we want to segment dark objects in a bright background since each 
object will be represented by a catchment basin. If we have objects with different intensities 
it is often better to apply the WT to the magnitude of the gradient instead Vu(x) since the 
gradient will be small inside each object. Figure 7.6a shows an image with several homogeneous 
patches. Figure 7.6b,b display the gradient magnitude and the watershed segmentation. 

How can we compute the watershed basins? we could simulate the water drop fall starting 
from each pixel of the input image. This was the first approach but it is very time consuming. 

A more interesting approach is based on immersion i.e., it assumes that the topographic 
surface is being flooded by water and the water level is increasing [9]. 

Let us assume that there is a hole located at each local minimum of the topographic surface. 
When the water level increases, the water starts to flood the deepest basins first and then all the 
others. Suppose the water in each basin has a different color. We want to build dams in such a 
way that the water from different basins does not mix and colors are preserved. The dams are 
built during the flood process. At the end of this process the dams define the watershed lines 
which define the segmentation of the image. Figure 7.7 illustrates the flood process. 

How can we implement dam construction in a computer? one way to do this in by precom- 
puting the set T p of all pixel locations x such that I(x) < p. This set defines all the regions 
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Figure 7.8: Region creation and merging 

which were flooded by water at level p. When we change from level p to level p + 1, two situa- 
tions must be detected (see Fig. 7.8): i) the presence of new connected components and ii) the 
merge of connected components (which should be prevented by building dams). 

The first situation can be easily detected by computing the connected components of T p+ i 
and checking if they all contain connected components of T p . If not a new basin is created. 
The second situation is also easily detected from the connected components of T n and T p+ \. If 
a connected component of T p+ i contained more than one connected component of T p we have 
a region merge which should be prevented. The dam can be constructed by performing a set 
of morphological dilations of the connected components of T n until they meet. In this process 
the connected components should not grow outside T p+ i. 

The constrained dilation of the regions C, C C T p is described by 

P+ V ) (7.21) 
C f ^T p+1 nD(C f ) 

where D(.) stands for the dilation operator. 



112 




Figure 7.9: Active contours 

7.7 Active Contours 

Active contours try a different approach instead of performing region labeling they try to fit 
the object boundary by a deformable curve. The curve is initialized in the vicinity of the object 
boundary and it is attracted toward the boundary by forces acting on the curve points (see Fig. 
7.9). 

The idea is very attractive. The main difficulty is how to attract the curve toward the object 
boundary (see Figure 7.10). How can this be done since we do not know the boundary points 
and we do not have a reliable criterion to do the match? 

To overcome this difficulty Kass et al. proposed the use of a potential function P(x) com- 
puted from the image (see Fig. 7. 10 (right)). P is chosen in such a way that intensity transitions 
become associates to valleys of the potential function. The contour model is then estimated by 
attracting each point toward the valleys. 

We will now formalize these ideas. Let x : [0, 1[— » M 2 be the parametrization of a curve in 
the image domain 1 . We associate an energy to each curve configuration 

E(x)= f 1 P(x(s))ds+ f 1 a\\x'(s)\\ 2 + f3\\x"(s)\\ 2 ds (7.22) 
Jo Jo 

where (.)' denotes a derivative with respect to s. The energy has two terms: an image dependent 
term (external energy) and regularization ter (internal energy) which tries to keep the curve 
derivatives small. This energy measures two conflicting goals. It measures the fitness of the 
model to the data and tries to keep the shape smooth. 

The minimization of E{x) is a problem of variational calculus. We wish to find the function 
x(s) which minimizes an integral of the form 

1 we assume in this section that the image domain is M? 
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Figure 7.10: Assignment problem: which boundary point should we choose? (left) potential 
function (right) 



Jo 



E= / F(x,x',x")ds 



(7.23) 



A necessary condition for optimality is given by the Euler equation. The Euler equation in 
this case is given by (see [7] and the proof in Appendix) 



ax"(s)-l3x""(s)-^-=0 
ax 



(7.24) 



This is a nonlinear differential equation whose solution is a stationary point of the energy. 
If we define the image force at the curve point x(s) as being the gradient of the potential at 
that point 



dP 

dx 



(7.25) 



and the internal forces as 



F int = ax"(s)-f3x""(s) (7.26) 
then the Euler equation defines an equilibrium of forces applied to each point of the curve 



Fimg(x) + F int (x) = 0 



(7.27) 



This should be true for the optimal contour configuration. 

To minimize E we can define a dynamic Euler equation which describes the evolution of the 
curve deformed by the internal and external forces. 

<9x 



dt 



^imgip^) ~l~ Ff n f;(x) 



(7.28) 



where x is a function of s and t. This equation usually converges to a stationary solution x(s) 
which verifies the (static) Euler equation as we wished. 
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To apply this equation we need to define the image potential. Several proposals have been 
made. Most of them try to make the potential small in the vicinity of the image transitions 
associated with the object boundary. Some choices are 



P(x) = -^(Mx)*u(x)) 

where h(x) is the impulse response of a low pass filter or 

P(x) = - J2 G CT (x-x 0 ) 

x 0 €T 

where T is the set of all the edge points in the image and 



G„(x) = 



27T<7 2 



e 2<r 



(7.29) 



(7.30) 



(7.31) 



is a Gaussian Kernel which spreads the influence of each edge point xq inside a circle. 



Discrete Snake 

The previous approach cannot be implemented directly since the curve is continuous and de- 
pends on an infinite number of parameters. To circumvent this difficulty the curve can be 
represented by a B-spline or approximated by a discrete sequence of points (see [6]). Let us 
consider the second case. The curve x : [0, 1[— » M 2 is replaced by a sequence of 2D points: 
x(0),x(l),...,x(iV-l). 

To estimate the curve configuration the dynamic Euler equation will be discretized. We will 
assume that t is an integer variable and replace the time derivative in (7.28) by a first order 
difference 

<9x(i, t) 



dt 



x(i, t + 1) - x(i,t) 



Therefore, 



x(i,£ + 1) = x(i,£) + F img (x(i,t)) + F int (x(i,£)). 



(7.32) 



(7.33) 



This equation is valid for all snake points i.e., for % — 0, . . . , N — 1. We can write this equation 
in a compact way using matrix notation. 

Let us define two matrices containing all the contour samples and the corresponding image 
forces, 

F im ,(x(0,t)) T 
Vim 9 (t)= : . (7.34) 



x(*) = 



x(0,t) T 
x(7V-l,t) T 



F im ,(x(7V_l,t))^ 
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x(t),Fi mg are 2 x N matrices. Then the snake update can be written as follows 

x(t + 1) = x(t) + Mx(t) + F imp (x(t)) 
where M is a circulant pentadiagonal matrix 2 



(7.35) 



M 



a 


b 


c 


0 .. 


. 0 


c 


b 


b 


a 


b 


c . . 


. 0 


0 


c 


c 


b 


a 


b .. 


. 0 


0 


0 



c 0 0 0 
b c 0 0 



b a b 
c b a 



(7.36) 



where a,b,c depend on the regularization parameters a,/3 as follows: a — —2a — 6/3, b — 
a + 4/3, c = — /?; M = 0 if no regularization forces are considered. 

Equation (7.35) provides a simple recursive rule to update the snake contour until it con- 
verges to a final configuration. The final configuration depends on the initialization. 

Figure 7.11 shows the estimation of the coronary vessels using snakes. These results were 
obtained using N — 30 samples and a — 1. Figure 7.11a shows the original contours. Figure 
7.11b displays the convergence of the snake algorithm and figure 7.11c displays the final contours 
when convergence is achieved. The final contour depends on the initialization as shown in this 
example. 

Snakes converges to the object boundary provided they are initialized close enough to the 
final configuration. They have however several limitations since they are attracted by other 
objects as well or by the clutter of the image background. They also perform poorly in the 
presence of concavities. These limitations can be partially overcome by using more sophisticated 
forces of improved shape models. For example if we have a statistical representation of the object 
shape derived from examples [8]. 



Discussion 

The snake converges to the object boundary if it is initialized close enough. However it has 
several limitations since it is also attracted toward other objects or to the background noise 
(clutter). These limitations can be aliviated by using better image forces and improving the 



2 a matrix is called circulant iif each line is obtained by a circular shift of the previous line. 
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Figure 7.11: Snake results: a) original image and initial contours; b) potential function and 
contour configuration at iteration 20, 50, 100, 150, 200 and c) final contour 

shape models. For example, statistical models of the desired shape can be used to make the 
shape estimate converge towards typical shape configurations [8]. 

Exercises 

• Derive an algorithm to update the discrete snake x = (x(0),x(l), . . . , x(7V — 1)) by in 
order to minimize the energy 

N-l N-l 

E=Y, Pi*®) + « E H x » - x (* - !)H 2 ( 7 - 37 ) 

i=0 i=0 

assuming that x(— 1) = x(iV — 1). Hint: the derivative of E with respect to x(p) must be 
zero. 

• Consider the discrete signal / = (4, 2, 1, 1, 3, 5, 4, 3, 3, 2, 2, 4, 5, 6). Compute the watershed 
transform suing the immersion algorithm. 

Appendix - Euler Equation 

The Euler equation is derived using variational calculus (e.g., see [7]). We will briefly derive 
the Euler equation for this case. 

Let x : [0,1 M 2 be a curve minimizing the energy (7.22). Let us modify x keeping the 
boundary conditions invariant. This can be done by adding a perturbation 

x(s) = x(s) +er](s) (7.38) 

where e is a constant and r\ is an arbitrary function such that 77(0) = 77(1) = 0, 77 (0) = 77 (1) = 0. 
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If e varies, the energy should have a minimum for e = 0 and 

= 0 (7.39) 



dE 
de 



6 = 0 

This stationary condition is the Euler equation. Let us now write it in a more explicit way 
taking the expression of E into account 

d_ r 1 

de , 



I' P(x + en) + £||x' + er/H 2 + |||x" + e V "\\ 2 ds = 0 
Jo ^ ^ 



(7.40) 



j 1 ^[P(x + erj) + |||x' + erj'f + |||x" + e V "\\ 2 } ds = 0 (7.41) 



Computing the derivative at e = 0 using the chain rule, 

T 



jf (£) ^ + a £V+^^V^ = 0 (7.42) 



A B 

Integrating by parts terms A and B we obtain 



f x ,t t/ ds = x ,T /y|J - x ,,T 7y d 5 (7.43) 
Jo Jo 

f x"V ds = x" T r,'\l - C x"'V ds (7.44) 
Jo Jo 

f 1 x" T rj" ds = x" T V '\l - x"' T V \l + f x"" T V ds (7.45) 
Jo Jo 



Replacing these integrals in (7.42) and using the initial conditions assumed we obtain 



C - ax" + W"] T n ds = 0 (7.46) 



to 

This equation is valid for any function r\. Therefore, (see [7]) 



d 4- - ax" + /?x"" = 0 (7.47) 
ax 



and the this the equilibrium equation used in the snake algorithm. 
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Chapter 



8 



Object 



Recognition 



8.1 Introduction 



Object recognition is a hot topic for more than three decades. Yest, it is not solved. Everyday 
thousands of images are produced in many contexts. Some of them are available on the net. 
We must have methods to access this information an efficient ways according to its context i.e., 
we need a visual google. Other recognition problems are signature recognition, handwritten 
letter recognition, faced recognition, car plate recognition and building recognition using aerial 
photographs. 

Some recognition tasks are simple. This happens if the object is easily detected in the image 
and can be easily characterized by its shape, motion or color. We have also extremely efficient 
algorithms to detect rigid objects with a known texture (e.g., building detais) invariant to the 
camera point of view and distance. However, it is still difficult to detect and recognize objects 
whose shape and texture presents a wide variety of different configurations (e.g., animals). Figure 
? shows several types of objects with different difficulty levels. 

Sometimes we do not want to recognize an object from a wide set classes. We just want to 
know if an object of a specific type (face, mine) is present in an image and where it is located. 
This problem is known as object detection. 

We can formalize the recognition and detection problems as follows: 

Problem 1: object recognition Given an image we and to recognize all the objects 
despite their position and orientation with respect to the camera. 
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Figure 8.1: Different object recognition challenges. 



A special case which occurs in many applications (e.g., surveillance and military applica- 
tions) is the detection of objects of interest (persons, planes, boats) in image and video. We 
can formulate that problem as follows. 

Problem 2: object detection Given an image we wish to detect the location of all the 
objects of a specific classe and their location. 

Object recognition problems may have a wide range of difficulties. The recognition of 2D 
objects with known shape and texture is easy even if the shape and texture information is 
randomly modified. There are also efficient techniques to deal with 3D objects with known 
texture observed from arbitrary points of view. Some of these techniques are rather recent 
[15]. The main challenge concerns 3D object recognition when the shape, color and texture 
dramatically change as it happens with animals (e.g., horses, dogs) and most objects in our 
daily life (cars, bycicles, flowers, chairs, tables, etc). 

Current systems are not able to solve the general recognition problem but they can provide 
acceptable solutions for many of the simpler problems. Let us discuss some examples. 

Examples 

License Plate Recognition 

Vehicles are identified by license plates. An automatic recognition of the license plate is a 
useful for surveillance and toll enforcement in highways. The image of the plate is obtained 
and processed to compensate for motion blurring and geometric deformation. Then charac- 
ter segmentation is performed and a character recognizer is used to retrieve the license plate 
number. 
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The last step is a typical Pattern Recognition problem which can be solved by standard 
techniques such as support vector machines, classification trees or neural networks. 

Detection of Moving objects 

It is often useful to detect moving objects (e.g., people, vehicles, boats) in images with cluttered 
backgrounds. This problem is solved by finding useful features (motion, color, texture) to 
discriminate the object from the background. 

One way to address this problem is by dividing the image into overlapping blocks and make 
a binary decision based on the block featires (motion, color, texture). This can be done by 
training a classifier algorithm as in the previous example. 

Silhuette Recognition 

The silhuette of an object (outer boundary) contains useful information to recognize the object. 
If the silhuette is a useful feature we can recognize the object using a 2D shape classifier. There 
are several shape features such as chain codes, invariant moments and Fourier descriptors, which 
can be used for this purpose. 

Face Detection 

Face detection is a very important problem in surveillance applications. The goal is to detect all 
faces present in an image. Current methods consider overlapping blocks and provide a binary 
decision for each block using a classifier. A sucessful example is [5]. 

Category Recognition 

This problem occurs when was wish to analyse large amounts of images from different sources 
(e.g., web). We wish to recognize the type of objects included in each image (e.g., vehicles, 
buildings, sky, animals). The main difficulty concerns the variety of shapes, texture and color 
within each class as well as the need for efficient methods able to cope huge amounts of infor- 
mation. 

8.2 Silluette Analysis 

The object silluette (boundary) is a useful cue for recognition. Several methods try to recognize 
objects from their silluette characterized by a closed curve 7 : [0, L[— » M 2 describing the object 
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boundary or by a binary image in which each point has a binary value 0, 1 depending on whether 
it belongs to the object or not. 

The same object may appear at different positions in the image and often with different 
scales and orientations. The recognition method should remain invariant or at least robust 
with respect to geometric transformations which may occur during the acquisition process. 

Region Comparison 

The simplest strategy consists of comparing two binary images directly on a pixel by pixel basis. 

Let / be a binary image representing an object of interest (/(x) = 1 iif x to the object) and 
let M be a binary image representing the object model. 

We can compare the images directly and check if there is a strong superposition between 
the regions marked with 1. Let A, B be the sets of image points with label 1 in the images /, 
M, respectively. The error between both sets of points can be measured by counting the pixels 
belonging to only one of the regions, and normalize this number by the object area 

#[(AuB)-(AnB)} 

D ~ iHAUW) (8 - 1} 

This method is simple but it has major drawbacks since it assumes that the objects are 
aligned. It does to deal with translations, rotations or scaling of the objects. 

Invariant Moments 

A binary object can be characterized by their central moments computed as follows. Let / be 
a binary image of the object. The center moment of order p, q is given by 

f JL PQ = ^2 (x ± - x 1 ) p (x 2 - x 2 ) q I(x u x 2 ) (8.2) 

( Xl ,x 2 )ez 2 

The center moments are invariant with respect to translation and they can be modified in order 
to enforce invariance with respect to rotation and scaling. 

Invariance with respect to scaling can be easily obtained by normalizing the center moments 
as follows 

- - ^ (8.3) 



,.(p+</)/2 
POO 



The invariance with respect to rotation is mutch harder but it can be achieved by the 



123 



nonlinear transformation proposed by Hu [1, 2] 

11 = 7720 + V02 

12 = (mo - V02) 2 + (2/yn) 2 

/3 = teo -3/7i 2 ) 2 + (3/721 -7703) 2 
14: = (7730 + V12) 2 + (mi + /fas) 2 

^5 = (7730 - 37712) (7730 + 7712) [(7730 + 7712) 2 - 3(7721 + 7703) 2 ] + (37721- (8.4) 
7703) (7721 + 7703) [3(7730 + 7712) 2 - (7721 + 7703) 2 ] 

16 = (7720 - 7702)1(7730 + 7712) 2 - (7721 + 7703) 2 + 47711(7730 + 7712) (7721 + 7703)] 

17 = (37721 - 7703X7730 + 77i 2 )[(t73o + 7712) 2 - 3(7721 + 7703) 2 ]- 

(7730 - 377i 2 )(772i + 77o3)[3(t73o + 7712) 2 - (7721 + 7703) 2 ]. 

This set of features are known as invariant moments and they have been used for object 
recognition. 

These features have an important drawback. Higher order models are very sensitive to small 
changes of the boundary since they lead to new terms multiplied by large factors (x — x) p (y — y) q . 
This effect amplifies small errors and makes the recognition more difficult. 

Chain Codes 

Uma representacao alternativa consiste em dividir o contorno do objecto 7 em trocos equiespacados 
e calcular a direccao de cada troco em relacao ao troco anterior. A mudanca de direccao e tipi- 
camente quantificada em 8 niveis 7ri/4, identificados pelo mdice mteiro i — 0, . . . , 7. A fronteira 
do objecto e assim caracterizada por uma sequencia de inteiros que representam as variacoes 
de direccao. Esta representacao e conhecida por codigo de cadeia diferencial (differential chain 
code) . 

O codigo de cadeia diferencial e invariante em relacao a translacoes e em relacao a rotacoes 
de multiplos de 7r/4. E ainda aproximadamente invariante em relacao a rotacoes de outras 
amplitudes. 

Uma dificuldade neste metodo consiste na escolha do ponto inicial a partir do qual se comeca 
a descrever o contorno. Essa escolha e por vezes feita calculando a interseccao do eixo de menor 
inercia do objecto com a fronteira e tomando um dos dois pontos de interseccao (p. ex., o mais 
afastado do centro de massa) como ponto inicial. 

O reconhecimento de objectos representados por codigos de cadeia passa a ser um problema 
de comparacao de sequencias que pode ser abordado por tecnicas variadas p. ex., alinhamento 
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temporal dinamico, cadeias de Markov nao observaveis ou metodos sintaticos de reconhecimento 
de padroes (gramaticas) . 

Descritores de Fourier 

O objecto pode ser reconhecido atraves da analise da fronteira. Seja j(s) uma curva fechada 
descrevendo a fronteira do objecto e ip(s) uma curva representando um modelo. 0 parametro 
s define univocamente cada ponto da curva. Por conveniencia admitiremos que estamos a 
trabalhar no piano complexo e 7(s),^(s) G C 

Ha muitas formas differentes de parametrizar uma curva. Ha ambiguidade na escolha do 
ponto inicial que corresponde a s — 0 e differentes formas de estabelecer a correspondencia 
entre os pontos da curva e os valores de s. Uma possibilidade e escolher s de forma a percorrer 
a curva com velocidade constante e de forma a que cada volta corresponda a um incremento de 
27r. As funcoes r y(s),ip(s) sao assim funcoes periodicas e podem por isso ser decompostas em 
serie de Fourier 

oo oo 
k= — oo k= — oo 

sendo os coeficientes da serie calculados atraves de produtos internos entre o sinal e as funcoes 
de base complexas (exponenciais). A norma L2 entre as duas funcoes 7(3), ip(s) pode ser obtida 
a partir dos coeficientes (relacao de Parseval) 

00 

llT-V'll 2 ^ E \ c k~dk\ 2 (8.6) 

k= — OO 

Contudo esta e uma ma medida para o reconhecimento pois requer que os objectos estejam 
alinhados. 

Vale a pena perceber qual a influencia das transformacoes geometricas sobre os coeficientes 
de Fourier. A tabela 8.1 resume este efeito para cada uma das transformacoes mencionadas. 
A translacao so altera o coeficiente de ordem 0 (DC). A rotacao altera a fase dos coeficientes 
somando-lhes um termo aditivo linear. As amplitudes fleam invariantes. A mudanca de ponto 
inicial tambem so altera a fase. A mudanca de escala altera as amplitudes mas de uma forma 
facil de compensar atraves de uma normalizacao. 

Um conjunto de descritores invariantes a estas transformacoes podem ser obtidos usando o 
modulo dos c/~ normalizado 

p* = j!f k \ , » ( 8 - 7 ) 

em que N e o numero de coeficientes usados para normalizacao. O coeficiente DC foi eliminado 
e os coeficientes de mdice negativo tambem devido a simetria. 
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Translacao 


7(s) + d 


c& + dJ(fc) 


Rotacao 


exp(jO)j(s) 


e j0 c k 


Ponto inicial 


70 - so) 


e jsok c k 


Escalamento 


cry(s) 


ac k 



(8.8) 



Table 8.1: Coeficientes de Fourier apos transformacoes geometricas 

8.3 Ajuste de template 

Um metodo simples para reconhecimento/deteccao de objectos consiste em obter uma imagem 
do objecto de interesse (template) T e compara-la com blocos de iguais dimensoes da imagem 
observada /. A comparacao e feita com base na norma l^. Assim, calcula-se 

£(x 0 ) = £(J(x-x 0 )-T(x)) 2 (8.9) 

X 

em que x 0 G Z 2 e o deslocamento aplicado a template. O objecto e detectado numa posicao x 0 
se E tiver um mmimo local nesse ponto e o mmimo for inferior a um limiar pre-deflnido 

£(x 0 ) < A (8.10) 

Este metodo designa-se por metodo de ajuste de template. A energia E e calculada tipicamente 
de forma esparsa ao longo de uma grelha de valores de x 0 que corresponde a ir deslocando a 
janela na imagem. Esta operacao deve ser feita com sobreposicao entre posicoes consecutivas 
da janela e e por isso pesada do ponto de vista computacional. 

O metodo de ajuste de template pode ser rescrito usando notacao matricial. Designando 
por t,x vectores obtidos por concatenacao das colunas de T e do bloco da imagem B vem 

£=||x-t|| 2 (8.11) 

O deslocamento x 0 esta implicito na escolha do bloco B. 

8.4 Detecgao Probabilistica 

Seja B um bloco de imagem que pode conter ou nao um objecto extraido de uma coleccao de 
objectos pre-defmidos Oi,i — 1, . . . , M e seja x um vector contendo todos os pixels de B. O 
vector x tern uma dimensao muito elevada (da ordem dos milhares). 
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O reconhecimento de Oi pode ser formulado como um problema de decisao. Dadas as 
distribuicoes do vector de observacoes para cada urna das classes (objectos) p(x.\Oi), classifica- 
se x na classe mais provavel 

P(Oi|x) = kp(x\Oi)P(Oi) i = 0, . . . , M (8.12) 

em que P{Qi\x) e a probabilidade a posteriori da classe Oi, P(Oi) a probabilidade a priori e 
k — l/p(x). Adicionou-se uma classe adicional Oo que designa a hipotese nula (nao ha objecto 
em B). 

Um caso particular eo do problema de deteccao em que ha apenas duas hipotese: O ou O 
(presenca ou ausencia do objecto de interesse em B). Neste caso pretende-se apenas saber se 
(deteccao de objecto) 

P(0|x) > P(0|x) (8.13) 

P(x|0)P(0) > P(x|0)P(0) (8.14) 

Por vezes conhecemos apenas a distribuicao de x quando ha objecto no bloco B. Neste caso, 
o teste pode ser aproximado por 

p(x|0) > A (8.15) 

em que A = p(x\0)P(O)\P(O). 

Se x tiver uma distribuicao gaussiana quando o objecto esta presente, x ~ N(x,R), entao 
a deteccao consiste em verificar se a desigualdade e verdadeira 

Ce -J(x-x) T il- 1 (x-x) > X ( 8<16 ) 

ou seja 

d(x) < a (8.17) 

d(x, x) = (x - xfR- 1 (x - x) (8.18) 

e a distancia de Mahalanobis e a — — 21og(A/C). 

O vector de media e a matriz de covariancia podem ser estimados a partir de um conjunto 
de treino de blocos {x% i — 1, . . . , N} por 

1 N 

x = ^E xZ ( 8 - 19 ) 

i=l 
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Figure 8.2: Estimagao de subespago 
1 N 

R = ^E( xi - i )( xi -*) T ( 8 - 2 °) 

Infelizmente as coisas nao sao tao simples. A dimensao muito elevada de x (~ 10 4 ) 
torna impraticavel o calculo da matriz de covariancia, nem e possivel reunir imagens de treino 
suficientes (>> 10 8 ) de forma a garantir uma boa estimativa de R. Estas dificuldades sao 
tratadas na seccao seguinte. 

8.5 Reconhecimento com PCA 

A analise em componentes principals e usada para reduzir a dimensao dos dados. 

Consideremos o caso mais simples em que se pretende aproximar um vector aleatorio xGl n , 
de media nula, por um vector pertencente a um subespaco com um unico vector de base, 
x = a(ft, </>Gl n (ver figura 8.2). Admitiremos que o vector (ft esta normalizado ou seja ||^|| = 1. 
O erro de aproximacao e e = x — acft. 

Pretende-se estimar a, (ft de forma a que o erro quadratico medio 

E = E{e T e} = E{x T x - 2o^ T x + a 2 (ft T (ft} (8.21) 

seja mmimo. O valor optimo para o coeflciente e a — x. T (ft e como (ft T (ft — 1 vem 

E = £{x T x - (ft T xx T (ft} = tr(R) - 0 T R(/> (8.22) 

em que R = i?{xx T } e a matriz de covariancia de x. A minimizacao de (8.22) sujeita a restricao 
de (ft T (ft = 1 obtem-se minimizando a funcao lagrangiana 

L = tr(R) - (ft T TL(ft + X((ft T (ft - 1) (8.23) 
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Calculando a derivada de L em ordem a (j> e igualando a zero obtem-se 

R0 = X(t> (8.24) 

que e uma equacao verificada pelos valores e vectores proprios da matriz R. Assim, cj> deve ser 
um vector proprio da matriz de covariancia R. Resta saber qual dos valores proprios devemos 
escolher para obter um mmimo. Substituindo (j> em (8.22) obtem-se 

E = tr(R) - A (8.25) 

Assim, a energia E e minimizada se (j) for o vector proprio associado ao maior valor proprio de 
R. Este vector corresponde a direccao de maior variacao dexee tambem conhecido por eixo 
de menor inertia. 

Consideremos agora o caso geral. Se quisermos aproximar x por um vector pertencente a 
subespaco de dimensao d < n o procedimento e identico ao anterior. O subespaco optimo e o 
subespaco gerado por d vectores proprios da matriz R associados aos maiores valores proprios. 
Assim 

d 

x = ^?/;^ (8.26) 

i=i 

com yi — x T 0f. Os vectores fa representam os principals modos de variacao de x e sao 
conhecidos por imagens proprias [4]. A equacao (8.26) traduz a decomposicao de x como uma 
combinacao de imagens proprias. 

E facil provar que o erro de aproximacao e uma combinacao linear dos restantes vectores 
proprios 

n 

e= Yl yrfi ( 8 - 27 ) 

i=d+l 

e a energia do erro e a soma dos n — d valores proprios mais pequenos ou seja, 

n 

E= Yj X i ( 8 ' 28 ) 

i=d+l 

admitindo que os valores proprios de R se encontram ordenados por ordem decrescente. 

Os coeficientes yi designam-se por componentes principals de x e a representacao de x como 
combinacao linear dos vectores proprios da matrix de covariancia designa-se por analise em 
componentes principals. 

A analise em componentes principals pode ser expressa em notacao matricial. Seja <l> G R nxn 

* = [^i...^n] (8.29) 
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a matriz de vectores proprios de R. Como os vectores proprios de uma matriz simetrica sao 
ortogonais (e foram normalizados) , a matriz <l> e unit aria ou seja, 

= I (8.30) 

e <l> -1 = <I> T . Por isso, a transformacao entre a imagem e os coeficientes pode ser escrita na 
forma 

x = *y y = $ T x (8.31) 

Esta transformacao e conhecida por transformacao de Karhunen-Loeve ou por analise em com- 
ponentes principais (PC A) [3]. 

A transformada de Kahunen-Loeve goza das seguintes propriedades. 

Transformacao unitaria: 

A matriz de transformacao <l> e uma matriz unitaria <l> _1 = <I> T . Por isso a transformacao 
diz-se unitaria. 

Invariancia da norma: 

A transformacao <l> preserva o comprimento dos vectores. Se y = <l>x, entao ||y|| = ||x||, Vx 
(||.|| designa a norma euclidiana). 
Diagonalizacao de R: 

A matriz de covariancia veriflca a propriedade 

* T R$ = A (8.32) 
em que A e a matriz diagonal definida pelos valores proprios da matriz R 

A = diag(X u . . . , A n ) Ai > A 2 > • • • > A n (8.33) 
Estimacao de subespaco: 

Pretende-se aproximar os vectores x G M n com covariancia R, por vectores (projeccoes) x 
pertencentes a um subespaco S C M n de dimensao d < n. O subespaco que minimiza o erro 
quadratico medio 

£ = £{||x-x|| 2 } (8.34) 

e o subespaco gerado gerado pelos vectores proprios associados aos m maiores valores proprios 
de R ou seja \I> = [3>i . . . &d] e 

x = *y y = ^ T x y eR d (8.35) 
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De entre todas as transformacoes ^ = [4>i . . . 4>d] e a transformacao optima no sentido de 
minimizar o erro quadratico medio de reconstrucao. Esta e uma importante propriedade da 
analise em componentes principais. As componentes principais sao as melhores variaveis para 
reconstruir o sinal a partir de informacao incompleta. 

O sinal x e decomposto em duas componentes ortogonais: sinal de aproximacao e erro 

d n 

x = x + e = ^ 2 / i * i + Yl V&i ( 8 ' 36 ) 

i=l i=d+l 

Como as duas componentes sao ortogonais, a energia de x e separada em duas partes (prove!) 

d n 

E{\\x\?} = E{\\x\\ 2 } + E{\\e\\ 2 } = £ Ai + E A * ( 8 - 37 ) 

1=1 i=d+l 

A energia de x e a soma dos primeiros m valores proprios e a energia do erro e a soma dos 
restantes 

n 

E{||x-x||}= X i ( 8 - 38 ) 

i=d+l 

O que fizemos ate ao momento foi decompor o espaco W 1 em dois subespacos ortogonais. 
Estes subespacos sao muitas vezes interpretados como um subespaco de sinal e um subespaco de 
ruido. O primeiro e gerado pelos d primeiros vectores proprios de R e o segundo pelos restantes 
vectores proprios. A dimensao d do espaco de sinal deve ser calculada de forma a que a energia 
do erro seja pequena (inferior a um limiar pre-flxado). 

Como veremos e facil estimar bem o subespaco de sinal mas e dificil estimar o subespaco de 
ruido. Isso tern consequencias como veremos. 

Vimos na seccao anterior que a deteccao de objectos e feita com base na distancia de 
Mahalanobis 

d(x) = j^R- 1 ^ = x T $A _1 $ T x = y T A _1 y (8.39) 

Assim, 

n 1 

d(x) = ^-y? (8.40) 

i=i 1 

A distancia de Mahalanobis tambem pode ser decomposta em duas partes correspondentes aos 
subespacos de sinal e de ruido 

<*(*) = Ex;»?+ E Y y " ( 8 - 41 ) 

i=l 1 i=m+l 1 

Em geral e facil calcular os valores proprios e os coeflcientes yi associados ao espaco de sinal 
(1° termo) mas nao ao espaco de ruido (2° termo). 
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Para ultrapassar est a dificuldade propoe-se em [4] que a distancia de Mahalanobis seja 
aproximada pela expressao seguinte em que o segundo termo e um valor aproximado do correcto 

rf w = E^ + i E yf ( 8 - 42 ) 

i=l 1 i=d+l 

em que A e a media aritmetica dos ultimos n — m — 1 valores proprios. A energia do erro pode 
ser facilmente obtida sem calcular os coeflcientes yi pois \\x — x|| 2 = ||x|| 2 — ||x|| 2 . Assim, 

d(x)=^ly? + i(|N| 2 -p|| 2 ) (8.43) 

1—1 1 

Este criterio de deteccao foi extensivamente testado em operacoes de deteccao de caras na base 
de dados FERET e comparado em concurso com as melhores tecnicas existentes na altura tendo 
obtido um dos melhores resultados [4]. 

Aspectos computacionais 

Nao e facil estimar R directamente a partir de uma sequencia de vectores x z devido a elevada 
dimensao dos vectores (~ 10 3 ). O numero de vectores de treino e geralmente da ordem das 
dezenas ou das centenas o que signiflca que a estimativa de R tern caracteristica nula (a maioria 
dos valores proprios sao nulos) e nao e possivel realizar uma analise de valores proprios com 
uma matriz de dimensoes tao grandes. Est as dificuldades sao facilmente ultrapassadas uma vez 
que se pretende calcular apenas os primeiros m valores e vectores proprios com m « n. 

Seja X = [x 1 ...x^], N < n uma matriz cujas colunas sao os vectores de treino. Os 
primeiros N valores e vectores proprios de R podem ser calculados a partir de uma matriz 
auxilizar de dimensao baixa N x n: ^X T X. 

Proposicao: [calculo dos valores e vectores proprios] 

Sejam (A^, v^),i = 1, . . . , N, a sequencia de valores e vectores proprios da matriz -^X T X. 
Entao (Ai, XvJ, i — 1, . . . , N sao valores e vectores proprios de R. 

A demonstracao e simples. Baseia-se na decomposicao em valores singulares da matriz 
de dados X = UDV T em que U, V sao matrizes ortogonais e D uma matriz diagonal cujos 
elementos da diagonal si sao os valores simgulares da matriz R. A matriz de covariancia 
exprime-se facilmente em termos destas matrizes 

R = ^XX T = ^UDV T VDU T = ^UD 2 U T (8.44) 

Comparando com a decomposicao em valores e vectores proprios da matriz (8.32) conclui-se 
que os valores proprios de R sao A^ = s 2 e os vectores proprios sao as colunas de U. 
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Repetindo o mesmo raciocmio com a matriz auxiliar 

^X T X = ^VDU T UDV T = ^VD 2 V T (8.45) 

conclui-se que a matriz auxiliar (de baixa dimensao) tern os mesmos valores singulares e portanto 
os mesmos valores proprios. Os vectores poprios da matriz auxiliar sao differentes no entanto 
podem ser relacionados com os de . 

Se Xi , forem um par de valor e vector proprio da matriz auxiliar 

X T X v, = A iVi (8.46) 

Multiplicando os dois membros por X vem 

XX T (X Vi ) = Ai(Xvi) (8.47) 

R(Xvi) = Ai(Xvi) (8.48) 
Ou seja Xv^ e vector proprio de R, como quenamos mostrar. 



8.6 Modelos de constelagao 

Os metodos anteriores sao aplicaveis quando os objectos sao faceis de segmentar ou tern uma 
aparencia visual com pouca variabilidade. Nao permitem resolver problemas gerais de recon- 
hecimento de objectos p. ex., reconhecer vacas, ou cadeiras ou arvores em posicoes gerais. Uma 
vaca pode ter varias cores possiveis. Pode estar em varias posicoes e ter por isso uma forma 
muito variada. Alem disso pode nao ser facil de segmentar. 

O problema de reconhecimento generico de objectos ou de classes de objectos tern um grande 
interesse e ainda um tema em aberto mas foram dados passos importantes nos ultimos 5 anos 
com vista a sua resolucao. 

A tendencia e a de desenvolver modelos locais para partes do objecto (chifre da vaca, 
barbatana do tubarao, olho de pessoa) ligados atraves de um modelo global que descreve a 
geometria relativa das varias partes (ver Figura 8.3). Esse modelo e conhecido por modelo de 
constelagao e tern sido desenvolvido por Perona e outros [6, 7]. 
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Figure 8.3: Modelos de partes (extraido de [8]) 
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Part III 
Vision 



Chapter 9 

Camera Model 



9.1 Introduction 

The principles of perspective were discovered during the Italian Renaissance by Brunelleschi, da 
Vinci and other architects and painters. They made drawings and paintings with a remarkable 
accuracy levels. They found how we see the 3D world. A famous work is Brunelleschi drawing 
of the church of Santo Spririto in 1436 (see Figure 9.1). 

Cameras project 3D points into the image plane. This transformation is not invertible 
unless we have additional information about the point location in space (e.g., depth). However, 
it allows us to obtain mutch information about a scene. Using a single camera it is possible to 
measure the length of a line on the floor, to determine the velocity of a car in a highway, to 
locate a football player on the field or to validate a goal line 1 . This section is inspired in two 
excellent textbooks [1, 2]. 

This chapter shows that camera performs a projection of 3-points in space onto the image 
plane described by a perspective projection 

x = nx (9.1) 

where X, x are the homogeneous coordinates of the points in space and in the image plane and 

n is a 3 x 4 matrix denoted as camera matrix. The camera matrix can be estimated from the 

data and fully characterizes the geometric transformation performed by the camera. 

1 recently, researchers from the university of Oxford showed that the third goal of England in the final of 
World Cup 1966 was not valid: the ball did not cross the line [3]. 
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Figure 9.1: Italian renaissance: Brunelleschi drawing (1436) and photo of the church of Santo 
Spirito 

9.2 Perspective Projection 

Digital cameras concentrate electromagnetic radiation (light rays) on the surface of the CCD 
sensor, using optical lenses. This process can be approximated by the pin hole model which 
assumes that the lenses are reduced to a small hole (optical center O) and all the light rays 
arriving at the sensor pass through that hole (see 9.2). Each point P in 3D space is projected 
into a point p obtained by intersecting the PO straight line with the image plane. 

The optical axis of the camera is a straight line orthogonal to the image plane which contains 
the optical center. The intersection of this line with the image plane is called the principal point 
of the image. 

Let us now consider the frames shown in figure 9.2: the XYZ frame centered at the optical 
center of the camera, and the image frame xy centered at the principal point of the camera. Let 
(X, Y, Z) be the coordinates of the point P in the camera frame and let (x, y) be the coordinates 
of the projected point p on the image plane. It can be easily shown by Tales theorem that 

x = -/§ y = -i\ (9.2) 

where / is the focal length of the camera. All these variables are measured in metric coordinates 
(m). This transformation which maps 3D points into image point is called the perspective 
projection. 

The images of the objects obtained using the pin hole model are inverted. To overcome this 
difficulty we can invert the image. We can do both transformation in a single step by assuming 
that the image plane is located in front of the optical center O (see Figure 9.3). This is called 
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Figure 9.2: Perspective Projection 



the frontal model. In this case, the coordinates of the projected point are given by 



X = fj 



y = f\ 



(9.3) 



This model is sometimes simplified by assuming / = 1 and in this case we call it normalized 
moid 

X Y , n t . 

x=j y=z (9-4) 

The perspective projection is a non linear transformation from the 3D space onto the image 
plane since the X, Y coordinates are divided by the depth. These equations become simpler if 
we represent the two points P,p using homogeneous coordinates. 

The homogeneous coordinates of two points P,p are defined by 



X = a 



X 
Y 
Z 
1 



x = /? 



X 

y 
l 



(9.5) 



where a, (3 are arbitrary constants. The homogeneous coordinates are obtained multiplying the 
Cartesian coordinates by a scale factor and adding the scale factor as an additional coordinate. 
The scale factor is arbitrary chosen. Therefore, two vectors of homogeneous coordinates X, X' 
represent the same point iif there is a scalar A / 0 such that X = AX ; . In this section we will 
use the same notation to denote a point written in Cartesian and in homogeneous coordinates. 
The meaning should be clear from the context. 

The normalized camera model in homogeneous coordinates is given by (check these equa- 
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tions) 



Figure 9.3: Perspective projection (front) 



X 
Y 
Z 
1 





X 




l 


0 
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0 



If we define a normalized projection matrix by 



n 0 = 



10 0 0 
0 10 0 
0 0 10 



= [io]e 



,3x4 



the camera model is given by 



Ax = n 0 x 

where A is a multiplicative constant (depth) . 



(9.6) 



(9.7) 



(9.8) 



9.3 Extrinsic Parameters 

Until now, we have chosen a 3D frame attached to the camera and located at the optical center 
of the camera (Figure 9.4). This hypothesis is not appropriate in practice. We can not easily 
measure the coordinates of point s with respect to the optical center hidden inside the camera. 

If we choose another frame with an arbitrary location and orientation, the new coordinates 
X 0 , are related to the camera coordinates by a rigid body transformation (rotation+translation) 
given by 2 

X = RX 0 + T (9.9) 



we are using Cartesian coordinates in this equation 
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Xot 



Figure 9.4: Extrinsic parameters 



where X 0 G M 3 is the vector of coordinated of the point measured in the world frame, R) is a 
rotation matrix, T G I 3 is a vector of translation. The parameters of the transformation (R, T) 
are called the extrinsic parameters of the camera since they depend on the camera location and 
orientation with respect to the world frame. The rotation matrix belongs to the special group 
of orthogonal matrices 50(3) i.e., it is an orthogonal matrix (R T R = I) and it preserves the 
orientation of the axis (det(R) = 1). 

The rigid body transformation can easily be written in homogeneous coordinates as follows 



X 
Y 
Z 
1 



R 
0 



T 

1 



X 0 
Y 0 

1 



(9.10) 



That is, 



X = gX 0 



g = 



R T 
0 1 

The camera model with extrinsic parameters is given by 



(9.11) 



Ax = IlogXo 



Ax = [R T]X 0 



(9.12) 



Example 

Consider a camera located in a rectangular room at the height h as shown in the figure. 
The optical axis of the camera belongs to the plane orthogonal to the wall and it make an angle 
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of —45° with respect to the horizontal plane. Given the world frame defined in the figure, we 
wish to compute the extrinsic parameters of the camera. 

Let us now solve this problem. Point p can be represented with respect to two frames: the 
room frame and the camera frame. These two sets of coordinates are related by X = RX 0 + T 
where (R, T) are the parameters of a rigid body transformation. 

The versors of the room frame can be expressed in the camera frame as follows 



ex 0 = e x 



e Yo — cos Oey 



suit 



e Zo — sin Oey — cos Oez 



(9.13) 



Therefore, 

1 0 0 

R= 0 cos 0 - sin 0 (9.14) 
0 sin 0 cos 0 

The translation vector should be chosen in such a way that the camera optical center has 
coordinates X 0 = [0, 0, h] T . Since RX 0 + T = 0, we obtain T = -h[0 sin 0 - cos 0] T . 



9.4 Intrinsic Parameters 

The metric coordinates (m) used so far to specify the location of point on the sensor surface 
are not practical. Image points are usually measured in pixels and frame is not centered at 
principal point of the camera. 

Therefore, the measured coordinates (V,?/) (pixels) are related to (x,y) (m) by 



x — fs x x -\- o x 



y' = fs v y + o y 



(9.15) 



where / is the focal length, s x ,s y are the conversion factors from meters to pixels and o x ,o y 
are the coordinates of the principal point in the new frame 3 . 

This transformation can be expressed in homogeneous coordinates as follows 



(9.16) 



x' 




fs x 0 o x 




X 


y' 




0 fSy Oy 




y 


l 




0 0 1 




l 



3 the focal length / is included in the transformation for the sake of convenience. 



142 



In practice we consider a more general transformation assuming that the entry 12 of the 
matrix may be different from zero. The transformation matrix becomes an arbitrary upper 
triangular matrix 

fs x fse o x 

x' = Kx K = 



f s v 



(9.17) 



0 

0 0 1 

The new parameter sq has also a geometric meaning: it is a diagonal scale (skew) factor. 

The parameters (fs x , fs y , fse, o x ,o y ) are denoted as intrinsic parameters of the camera since 
they only depend of the camera and are independent of the camera position and orientation. 
Matrix K G l 3x3 is called the matrix of intrinsic parameters 4 . 

The camera model with intrinsic and extrinsic parameters is given by 



Ax' = KllogXo 
Ax' =K[R T]X 0 

The camera model can also be expressed as 

Ax' = nx 0 



(9.18) 
(9.19) 

(9.20) 



where nEM 3x4 isa3x4 matrix denoted as the camera matrix. For the sake of simplicity we 
will replace Xq by X from now on. 



9.5 Camera Model in Cartesian Coordinates 



Let us consider the camera model (9.20). The camera matrix n can be written as 



n = 



T 
^2 



where irf is the i-th line of n. Using (9.20) we obtain A = ttJX 0 and 



7rfX 0 



y 



ttJXq 



(9.21) 



(9.22) 



7T^X 0 7T§ 

this is a nonlinear transformation from the 3D space to the image plane known as the perspective 
projection. 



4 a zoom operation changes matrix K. 
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Figure 9.5: Projection of a 3D line: a) general case, b) parallel to image plane; c) containing 
the optical center 

The perspective projection was discovered by the Greeks and later in Italy by the renaissance 
architects Brunelleschi and Leonardo da Vinci. It allows an accurate account of transformation 
produced by the eye when we observe a scene and wish to draw it. 

9.6 Points, lines and planes 

The projection of a point in space P onto the image plane is a point p whose coordinates are 
given by (9.20) (see Fig. 9.3). Let us now consider the projection of lines and planes. These 
projections were heavily used in painting and architecture since the Italian Renaissance with 
the works of Brunellesci and Leonardo da Vinci 

The projection of a 3D line r on the image plane is a line or a point (see Fig. 9.5). The 
projection is a line if the line r does not contain the optical center (Fig. 9.5a,b). If the 3D line 
contains the optical center the projection is a point (Fig. 9.5c). 

Let us discuss the first case and let us consider a point moving along the straight line. The 
projection of the point on the camera moves along the projected line and converges to a point 
called vanishing point. The vanishing point can is obtained by intersecting the image plane 
with a line parallel to r and containing the optical center O (see Figure 9.5a). 

How can we compute the coordinates of the vanishing point? Let us consider a straight line 
r in the 3D space. The straight line is the set of all points 

X = X 0 + aT (9.23) 

where X = [X Y Z 1] T is the position of a point on the straight line in homogeneous coordinates, 
X 0 = [Xq Y 0 Z 0 1] t is the initial point, T = [T x T y T z 0] T is the tangent vector and a is a free 
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Figure 9.6: Projection of a cubic house using 3 vanish points. 



parameter (time). Note that the fourth coordinate of T is zero. The projection of X on the 
image plane is (see (9.20)) 

x = nX 0 + aUT (9.24) 

The point x belongs to the straight line defined by the points xi = nX 0 ,x 2 = IIT (in 
homogeneous coordinates). When a tends to infinity 

x = IIT (9.25) 

The vanishing point in Cartesian coordinates is given by 

x = 7j7~~~ y = -\- (9.26) 

with T = [T x T y T z 0] T . Several conclusions can be draw. The vanishing point does not depend 
on the position of the line in space (X 0 ). It only depends on the orientation. Therefore a set 
of parallel lines has the same vanishing point. This is a well known property of the perspective 
projection. There is one vanishing point for each set of parallel lines. For example, if we want 
to draw a building with cubic shape we must consider three vanishing points. Furthermore, the 
vanishing point of all lines belonging to a plane are located on an image line called the horizon 
(Fig. 9.6). 

The projection of a plane is a plane except if the plane contains the projection center. In 
this case the projection is a straight line. Let us briefly discuss the horizon line we observe 
when we take a photograph of the sea. Is the horizon caused by the curvature of the earth 
(see Fig. 9.7)? Try to convince yourself that you would see a line of horizon even the case was 
plane. 
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O I horizon? 




Figure 9.7: Is the line of horizon produced by the curvature of earth. 

9.7 Camera Calibration 

An important issue concerns the estimation of the camera parameters from an image. This 
process is known as the camera calibration. 

Camera calibration is often done using a calibration object with with known geometry (see 
figure 9.8). The calibration object should contain a set of easily detected points. We then 
obtain an image of the calibration object and accurately detect all the calibration points. The 
calibration data is a set of 3D points and their projections 

T={(X i ,x i ),i = l,...,n} (9.27) 

where X J Gl 3 ,x J E M 2 are the coordinates of the i-th point in space and in the image plane. 

The camera parameters can be chosen in order to minimize the distance from the points 
projected by the camera and by the model. This can be done minimizing the quadratic cost 
functional 

n 

^^^-/(XV))! 2 (9.28) 

where / : M 3 — > R 2 is the perspective projection (9.29) and tt are the model parameters. The 
minimization of (9.28) is a non convex optimization problem which can be numerically solved. 
However, the cost function has many local minima and the optimization algorithm can be 
trapped by one of them. 

To overcome this difficulty simpler methods are used based on a different optimization 
criterion. The camera model (9.29) can be written as follows 

XX T 7T 3 = X T 7Ti 2/X T 7T 3 = X T 7T 2 (9.29) 

where 7^ is the i-th line of matrix II. 
Using matrix notation we obtain, 

" X T 
0 



0 



-xX T 



7T = 0 



(9.30) 
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Figure 9.8: Calibration object 



where tv G M 9x1 is a vector with the coefficients of matrix <l>. 

Each pair of points (X, x) leads to 2 linear equations. Therefore, the data set T defines a 
system of 2n equations 

Mtt = 0 



(9.31) 



where 



M: 



X 1T 
0 



0 

X 1T 



-y^ 1T 



0 



0 



2nx9 



(9.32) 



-x n X nT 

_ y n X nT 

We conclude from (9.31) that tv should belong to the null space of M. In practice M has rank 
9 and tv is obtained solving (9.31) by least squares. This is done minimizing the norm of the 
residue 

e = Mtt (9.33) 

assuming that the camera coefficients are normalized i.e., ||7r|| = 1. This is an optimization 
problem with constraints. We can solve this problem using a Lagrangean function 



L = 7r T M T M7r + A(e T e - 1) 



(9.34) 



The minimization of L can be done computing the derivative and making it equal to zero. We 
obtain 

M t Mtt = -Att (9.35) 

Therefore tv must be an eigen vector of M T M. It can be easily shown that tv is the eigenvec- 
tor associated to the smallest eigen value of M T M. The estimation of tv is done by simply 
computing the eigendecomposition of matrix M T M G M 9x9 . 



147 



After obtaining II we still need to obtain the intrinsic and extrinsic parameters of the 
camera. This is possible and can easily be done. Let II = [A b] where A is a square matrix 
formed by the first three columns of II and b is the las column. Therefore, 

A = KR b = KT (9.36) 

where K is an upper triangular matrix and R is a rotation matrix. 

The QR decomposition from Linear Algebra solves a similar problem. It decomposes a 
square matrix M into the product 

M = QR (9.37) 

where Q is a rotation matrix and R an upper triangular matrix. We cannot use this result 
directly because the upper triangular and the rotation matrix appear in a different order. This 
difficulty is overcome by inverting A. Therefore, 

A" 1 = R^K 1 (9.38) 

where R _1 is a rotation matrix and T _1 is an upper triangular matrix 5 . We can apply the 
QR decomposition to A -1 to obtain R _1 , K -1 . The translation vector is then easily computed 
T = K _1 b. The intrinsic and extrinsic (pose) parameters of the camera are easily obtained in 
this way. 



5 the inverse of a rotation matrix is still a rotation matrix and the inverse of an upper triangular matrix is 
still upper triangular 
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Chapter 10 

Shape and Motion Estimation 



10.1 Introduction 

Cameras project the information of 3D world into 2D images. There is loss of information in 
this transformation (depth and occlusions). Can we invert the projection performed by the 
camera and recover 3D shape and motion? This is the question addressed in this chapter and 
it is fundamental issue if we want to use cameras to measure the world. 

The main question we will address is: can we recover the 3D geometry and motion from 
images produced by a moving camera? 

It is difficult to answer this question using a single image. The scene points are projected 
onto the image plane by a perspective projection (Figure 9.3). The projection of a 3D point 
P defines an optical ray but does not convey any information about depth. Therefore, the 
projection is not enough to reconstruct P 1 . 

Suppose we have two views of the same scene with overlap (some points are visible in both 
views, see Figure 10.1). Knowing the projections pi,P2 of a point P in both views and knowing 
the optical centers of the cameras 0\ , O2 , it is possible to recover P computing the intersection 
of the lines 0\p\ 1 02P2 (see figure 10.2). This procedure is called triangulation. 

We can repeat this procedure for many 3D points and estimate the objects surfaces by 
interpolation. This approach is rather general since it make no assumptions about the object 
shape. However it has some limitations. It requires accurate models of the cameras and assumes 
that we can detect many pair of corresponding points in both cameras. These restrictions can 
1 in some cases, depth can be recovered from a single image (e.g., if P belongs to a known plane) 
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Figure 10.1: Pair of images 
P 




Figure 10.2: Triangulation: point P is obtained by intersecting the lines P\0\ 1 p202 

be aleviated using dense approach but we will not pursue that line here. 

The estimation of a 3D point requires information about the intrinsic parameters of the 
camera and the camera motion. Three cases can be considered: i) known internal model and 
camera motion; ii) known internal model and unknown camera motion; iii) unknown internal 
model and camera motion. Table 10.1 summarises these three cases. 

The first two problems are known as 3D reconstruction with calibrated camera. The last 
problem is more difficult and it is known as the uncalibrated reconstruction. 

10.2 Epipolar Geometry 

The reconstruction of P is obtained by intersecting the corresponding optical rays (see Figure 
10.2). However, they may not intersect in space. This happens very often since the image 
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Available information 


Information to be estimated 


Intrinsic parameters 
motion 


shape 


intrinsic parameters 


shape and motion 




shape and motion 



Table 10.1: Hierarchy of reconstruction problems 

P 




Figure 10.3: Epipolar geometry 

coordinates are corrupted by noise and the camera model is approximate. This raises a question: 
how can we guarantee that two optical rays associated to a pair of views intersect in space? 

The answer to this question is based on the concept of epipolar geometry and involves three 
concepts: epipolar plane, epipolar line and epipoles. 

The optical centers 0\ , O2 , p and the scene point P define a plane in 3D space known as the 
epipolar plane associated to P. The optical rays belong to the epipolar plane. Furthermore, 
the epipolar plane intersects the image planes in two lines h^h known as epipolar lines. 
Therefore, each point P in space defines a pair of epipolar lines. These lines have an important 
role: two points pi,P2 are projections of a point P if and only if they belong to corresponding 
epipolar lines. Therefore, if we search for a pair of corresponding points we should restrict the 
search to pairs of epipolar lines. This allows important computational savings. 

There are two points which belong to all the epipolar lines. These points are known as 
epipoles. The epipoles ei, e2 are obtained by intersecting the line O1O2 with the image planes. 
There are degenerated cases in which the epipolar plane is not uniquely defined. If 0\ — O2 
(pure rotation) there are an infinite number of epipolar planes. 

We will now show that the epipolar geometry can be estimated from a pair of images without 
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knowledge about the camera motion or the 3D scene. 



10.3 Essential Matrix 

Let us for the sake of simplicity that the cameras are calibrated. This means that we know the 
intrinsic parameters of the camera. The projected points xi,x 2 can be pre-processed in order 
to compensate the intrinsic parameters of the camera. Pre-processing is done by multiplying 
the homogeneous coordinates of the points by the inverse of matrix K: x — > K _1 x. 
Under this hypothesis, the camera model in both positions is very simple 

Aixi = Xi A 2 x 2 = X 2 (10.1) 

where xi , x 2 G M 3 are the homogeneous coordinates and Xi , X 2 G M 3 are the Cartesian co- 
ordinates of P measured in the cameras frames. Since X 2 = RXi + T where (R, T) are the 
parameters of a rigid body transform, we obtain 

A 2 x 2 =RAiXi +T (10.2) 

This equation creates a constraint on the projected points xi, x 2 . It is not easy to interpret 
this equation since it depends on the scale factors Ai,A 2 which are unknown. We will modify 
the equation to get rid of the scale factors. 

Let us compute the external product of both members of (10.2) with T = [tit 2 ts] T . This 
can be done by multiplying both members by the skew symmetric matrix 



T = 



0 -t 3 t 2 

h o -h 
-t 2 h o 



(10.3) 



We obtain 

A 2 f x 2 = Aif R Xl (10.4) 

since T x T = 0. If we compute the inner product of both members of the equation by x 2 we 
obtain 

x^f Rxi = 0 (10.5) 

since Tx 2 is orthogonal to x 2 . This is a constraint that each pair of projected points (xi,x 2 ) 
must obey. This is also a sufficient condition for guarantying that two points in the image 
planes of both cameras are the projections of the same 3D point P. 
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The matrix E = TR is called the essential matrix and plays an important role in 3D 
reconstruction from pairs of images. Furthermore, the matrix E can be directly estimated from 
a pair of images and allows us to compute epipolar geometry: epipolar lines and epipoles. from 
a pair of images and allows us to compute epipolar geometry: epipolar lines and epipoles. 

Before we do this we need to see how to define a straight line in homogeneous coordinates. 
A straight line is defined by the equation 

l T x = 0 (10.6) 

where 1 E M 3 is a vector containing the unit norm and the distance to the origin. The epipolar 
lines can now be easily computed from the essential matrix (10.5). If we keep xi (X2) constant 
we obtain the equation of a line for the variable x 2 (xi). This is the epipolar line. Both epipolar 
lines are defined by 

If xi = 0 1^x2 =0 h = E T x 2 1 2 = Exi (10.7) 

The epipoles can be easily computed as well. They must belong to all epipolar lines. There- 
fore, 

Eei = 0 E T e 2 = 0 (10.8) 

The epipole ei is associated to the null space of the essential matrix E and the epipole e 2 is 
associated to the null space of E T . Each of them can be easily computed from E. If we use the 
equation E = TR we can also conclude that ei = R T T e e 2 = T (apart from a scaling factor). 

Matrix E = TR has interesting properties. It can be shown that E is an essential matrix 
iif it has a null eigen value and the other two are equal [1]. The eigen decomposition of E is 
given by 

E = UAV T A = diag{a, a, 0} (10.9) 

The space of matrices which verify this equation is called the essential space. 

We will now consider an important question: is it possible to invert the application (R, T) — > 
E = TR. The answer is yes apart from the fact that there are two possible solutions for the 
inverse transform. Let E be a matrix belonging to the essential space, then E = TR with: 

R x = UR Z (+ J)EU T Ti = UR Z (+ J)V T (10.10) 
R 2 = UR Z (-|)SU T T 2 = UR Z (-|)V T (10.11) 
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where R^(^) is a rotation of amplitude 9 with rotation axis z 

R*(0) = 



cos(<9) - sin(0) 0 
sin(<9) cos(<9) 0 
0 0 1 



(10.12) 



This result has an important impact in practical applications. Motion information can be 
recovered from the essential matrix! 



10.4 Estimation of the Essential Matrix 

How is the essential matrix estimated? 

The essential matrix can be estimated from a pair of images by linear estimation algorithms. 
The fundamental equation (10.5) is linear in the parameters of E. Defining xi = [x\ , y\ , zi], x 2 = 
[x 2j y 2j z 2 ] we obtain 

[x ± x 2 x x y 2 xiz 2 y±x 2 y±y 2 y\z 2 z x x 2 z x y 2 Zl z 2 ]e = 0 (10.13) 

where e = [en e 2 ± . . . ess] T is a vector with the elements of the essential matrix. The previous 
equation can be written as follows 

a T e = 0 a = xi(g)x 2 (10.14) 

where 0 stands for Kronecker multiplication of two matrices 2 . 

Each pair of points defines a linear equation for e. If we know n pairs of points detected in 
both images {(x^,X2, )} we obtain a system of equations 

Me = 0 (10.15) 

where 

M G M nx9 (10.16) 

Vector e belongs to the null space of M. 

Since e is defined up to a scale factor we assume that ||e$ = 1. It is not possible to solve 

equation (10.16) in an exact way with the constraint. If we have more than 8 points, matrix 

2 M = A(g)Bisa block matrix M = [M^] where each block is the product of an element of A by B as 
follows NLn = aa'B. 



M = 
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M will have rank 9 and the null space has a single element 0).To overcome this difficulty we 
solve the problem using the least squares method by minimizing the sum of the residues 

£=||Me|| 2 (10.17) 

with the constraint ||e|| 2 = 1. This problem can be solved using the Lagrangean function 

E= ||Me|| 2 -A(||e|| 2 -l) (10.18) 

E = e T M T Me - A(e T e - 1) (10.19) 
Computing the derivative of E with respect to e we obtain 

M T Me - Ae = 0 => M T Me = Ae (10.20) 

this equation is verified if e is an eigenvector of M T M. The minimum of E is achieved if we 
choose the eigenvector associated to the smallest eigenvalue of M T M. 

In general, the least squares estimate is not a valid essential matrix since it does not belong 
to the essencial space. The next step consists of approximating the least squares estimate by 
a valid essential matrix. This is done projecting the estimate onto the essential space. Let the 
SVD decomposition of E be 

E = UAV T A = diag{al, a2, <t3} (10.21) 

The projection of E on the essential space is 

E' = UA'V T A' = diag{a, a, 0} (10.22) 

em que a — (oi + cr 2 )/2. 

This algorithm is known as the 8 point algorithm since it requires at least 8 pairs of points. 
This algorithms was proposed by Longuet-Higgins (Psichologist, Cambridge) in the scope of 
the study of human perception and it is summarized in Table 10.2. This algorithm is very 
sensitive to measurement errors but can be significantly improved if the data is preprocessed. 
The mass center of the image points should be translated to the origin and the variance should 
be normalized. 

10.5 Shape and Motion estimation (Calibrated camera) 

Motion estimation can be done using the essential matrix E using (10.10). We will now focus 
on the estimation of shape of the objects viewed by both cameras. We will assume that all the 
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8 point algorithm (Longuet-Higgins) 



1. Data acquisition: determine N pairs of corresponding points in both cameras: 
{(xi,xj),i = l,...,iV} (N>8) 

2. Pre-processing: center and normalize each set of points x^x^ subtracting their mean 
and scaling them by the inverse standard deviation. These operations will be denoted 

by Ai, A 2 . 



3. Least squares estimate: compute the eigenvector associated to the smallest 
eigenvalue of M T M. Matrix M is defined in (10.16). This estimated is denoted by E 

4. Projection: determine the SVD of E E = UAV T , 

A = diag{al, cr2, aS} and compute the projection of E the essencial subspace 
E' = UA'V T , em que A' = diag{a, cr, 0} e a — (01 + cr 2 )/2. 

5. Post-processing: E = A^E'Ai 

Table 10.2: 8 point algorithm for the estimation of E 
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3D points to be reconstructed are independent and will be estimated in an independent way 3 . 

There are several ways to estimate P from the projections xi,x 2 . They are all equivalent 
when the cameral model is perfect and the observation noise is zero. However, the methods 
perform in different ways when applied to real data. We will now discuss a simple linear 
algorithm to estimate the 3D position of observed points. 

If we adopt the model (10.1) we conclude that the estimation of P is equivalent to the 
estimation of Ai or A 2 . In general it i not possible to enforce both conditions in an exact way. 
We prefer instead to define algebraic errors t\ , e 2 G M 3 between both members of the equations 

Cl = Aixi - X e 2 = A 2 x 2 - (RX + T) (10.23) 

Since it is not possible to guarantee that all the errors are zero, the vector X is estimated by 
minimizing a quadratic cost functional 

E = ||AiXi - X|| 2 + ||A 2 x 2 - RX - T|| 2 (10.24) 

where ||.|| is the Euclidean norm. We will minimize 

E = (Aixi - X) T (AiX! - X) + (A 2 x 2 - RX - T) T (A 2 x 2 - RX - T) (10.25) 

A necessary condition verified by the minima of E is 

Computing the derivative of E with respect to X (see Appendix) 

-(Aixi - X) - R T (A 2 x 2 - RX - T) = 0 (10.27) 

(I + R T R)X = Aixi + A 2 R T x 2 - R T T (10.28) 
This is a system of 3 linear equations. The solution is given by 

X = (I + R T R)~ 1 (AiXi + A 2 R T x 2 - R T T) (10.29) 

The matrix inversion can be done once for all the data points since it does not depend on the 
points detected in the images. 

This method is algebraic and the cost functional does not have a geometric meaning. To 
overcome this difficulty we could reconstruct P by estimating the point closest to the straight 
lines oixi, o 2 x 2 . This method is also linear. 



3 some methods take the opposite approach and try to estimate a dense depth map 
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Let us consider a third approach based on the image plane. Let us assume that the projected 
pints are corrupted by noise. Therefore the available measurements are 

x' x = xi + wi x 2 = x 2 + w 2 (10.30) 

where wi , w 2 is the observation noise and xi , x 2 are the true projections which should belong to 
a pair of epipolar lines and should verify (10.5). The reconstruction problem can be formulated 
as follows O problema: what is the pair of epipolar lines closest to x^x^? the answer to this 
question involves the solution of a nonlinear optimization problem. 

Exercises 

1. [Reconstrucao 3D] Let £i,x 2 be the projections of a point P by a calibrated camera at 
two positions and let (R, T) be the motion parameters. Compute the coordinated of 
the space point X closest to the straight lines 01X1, o 2 x 2 . (Hint: use (10.1) and motion 
transformation) 
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Apendices 



Derivadas em ordem a matrizes 

A derivada de uma funcao escalar /(X) em ordem a uma matriz X de dimensoes n x m e uma 
matriz defmida por 

df ^ - 9f (10.31) 



dxy 4i dXij 

Nalguns textos a matriz de derivadas aparece transposta em relacao a esta definicao aqui adop- 
tada. 

Em seguida mencionam-se algumas propriedades do traco e do determinante de uma matriz 
admitindo-se que X, A e B sao matrizes de dimensoes apropriadas. 
Propriedades do traco: 

±tr{X} = I T (10.32) 

-^rfr{AXB} = A T B T (10.33) 

-^-ir{AX T B} = BA (10.34) 

-^-tr{XX T } = 2X (10.35) 

-^-tr{AXBX} = A T X T B T + B T X T A T (10.36) 

-^-tr{AXBX T } = AXB + A T XB T (10.37) 

-^7*r{AX -1 B} = -X _1 BAX _1 (10.38) 

d 'X| = |X| X" T (10.39) 



dX 
d_ 
dX 



d log|X|=X- T (10.40) 



-^r | AXB | = | AXB | X" T (10.41) 
^-\X\ n = n\X\ n X~ T (10.42) 
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